diff --git a/nomad/datamodel/metainfo/public.py b/nomad/datamodel/metainfo/public.py index cc22bfc99e7c8e61e4ac3bc90d26c0cbe1f1ac8e..b3a70193473c4dff837cd01e7995b53a1f720d3f 100644 --- a/nomad/datamodel/metainfo/public.py +++ b/nomad/datamodel/metainfo/public.py @@ -2,7 +2,7 @@ import numpy as np # pylint: disable=unused-import import typing # pylint: disable=unused-import from nomad.metainfo import ( # pylint: disable=unused-import MSection, MCategory, Category, Package, Quantity, Section, SubSection, SectionProxy, - Reference + Reference, MEnum ) from nomad.metainfo.legacy import LegacyDefinition @@ -2025,111 +2025,6 @@ class section_gaussian_basis_group(MSection): a_legacy=LegacyDefinition(name='number_of_gaussian_basis_group_exponents')) -class section_k_band_normalized(MSection): - ''' - This section stores information on a normalized $k$-band (electronic band structure) - evaluation along one-dimensional pathways in the $k$ (reciprocal) space given in - section_k_band_segment. Eigenvalues calculated at the actual $k$-mesh used for - energy_total evaluations, can be found in the section_eigenvalues section. - ''' - - m_def = Section(validate=False, a_legacy=LegacyDefinition(name='section_k_band_normalized')) - - k_band_path_normalized_is_standard = Quantity( - type=bool, - shape=[], - description=''' - If the normalized path is along the default path defined in W. Setyawan and S. - Curtarolo, [Comput. Mater. Sci. **49**, 299-312 - (2010)](http://dx.doi.org/10.1016/j.commatsci.2010.05.010). - ''', - a_legacy=LegacyDefinition(name='k_band_path_normalized_is_standard')) - - section_k_band_segment_normalized = SubSection( - sub_section=SectionProxy('section_k_band_segment_normalized'), - repeats=True, - a_legacy=LegacyDefinition(name='section_k_band_segment_normalized')) - - -class section_k_band_segment_normalized(MSection): - ''' - Section collecting the information on a normalized $k$-band segment. This section - stores band structures along a one-dimensional pathway in the $k$ (reciprocal) space. - - Eigenvalues calculated at the actual $k$-mesh used for energy_total evaluations are - defined in section_eigenvalues and the band structures are represented as third-order - tensors: one dimension for the spin channels, one for the sequence of $k$ points for - the segment (given in number_of_k_points_per_segment), and one for the sequence of - eigenvalues at a given $k$ point. The values of the $k$ points in each segment are - stored in band_k_points. The energies and occupation for each eigenstate, at each $k$ - point, segment, and spin channel are stored in band_energies and band_occupations, - respectively. The labels for the segment are specified in band_segm_labels. - ''' - - m_def = Section(validate=False, a_legacy=LegacyDefinition(name='section_k_band_segment_normalized')) - - band_energies_normalized = Quantity( - type=np.dtype(np.float64), - shape=['number_of_spin_channels', 'number_of_normalized_k_points_per_segment', 'number_of_normalized_band_segment_eigenvalues'], - unit='joule', - description=''' - $k$-dependent energies of the electronic band segment (electronic band structure) - with respect to the top of the valence band. This is a third-order tensor, with - one dimension used for the spin channels, one for the $k$ points for each segment, - and one for the eigenvalue sequence. - ''', - a_legacy=LegacyDefinition(name='band_energies_normalized')) - - band_k_points_normalized = Quantity( - type=np.dtype(np.float64), - shape=['number_of_normalized_k_points_per_segment', 3], - description=''' - Fractional coordinates of the $k$ points (in the basis of the reciprocal-lattice - vectors) for which the normalized electronic energies are given. - ''', - a_legacy=LegacyDefinition(name='band_k_points_normalized')) - - band_occupations_normalized = Quantity( - type=np.dtype(np.float64), - shape=['number_of_spin_channels', 'number_of_normalized_k_points_per_segment', 'number_of_normalized_band_segment_eigenvalues'], - description=''' - Occupation of the $k$-points along the normalized electronic band. The size of the - dimensions of this third-order tensor are the same as for the tensor in - band_energies. - ''', - a_legacy=LegacyDefinition(name='band_occupations_normalized')) - - band_segm_labels_normalized = Quantity( - type=str, - shape=[2], - description=''' - Start and end labels of the points in the segment (one-dimensional pathways) - sampled in the $k$-space, using the conventional symbols, e.g., Gamma, K, L. The - coordinates (fractional, in the reciprocal space) of the start and end points for - each segment are given in band_segm_start_end_normalized - ''', - a_legacy=LegacyDefinition(name='band_segm_labels_normalized')) - - band_segm_start_end_normalized = Quantity( - type=np.dtype(np.float64), - shape=[2, 3], - description=''' - Fractional coordinates of the start and end point (in the basis of the reciprocal - lattice vectors) of the segment sampled in the $k$ space. The conventional symbols - (e.g., Gamma, K, L) of the same points are given in band_segm_labels - ''', - a_legacy=LegacyDefinition(name='band_segm_start_end_normalized')) - - number_of_normalized_k_points_per_segment = Quantity( - type=int, - shape=[], - description=''' - Gives the number of $k$ points in the segment of the normalized band structure - (see section_k_band_segment_normalized). - ''', - a_legacy=LegacyDefinition(name='number_of_normalized_k_points_per_segment')) - - class section_k_band_segment(MSection): ''' Section collecting the information on a $k$-band or $q$-band segment. This section @@ -2211,6 +2106,63 @@ class section_k_band_segment(MSection): a_legacy=LegacyDefinition(name='number_of_k_points_per_segment')) +class section_band_gap(MSection): + ''' + This section stores information for a band gap within a band structure. + ''' + m_def = Section(validate=False, a_legacy=LegacyDefinition(name='section_band_gap')) + + value = Quantity( + type=float, + unit="joule", + description=""" + Band gap energy. + """, + a_legacy=LegacyDefinition(name='value') + ) + type = Quantity( + type=MEnum("direct", "indirect"), + description=""" + Type of band gap. + """, + a_legacy=LegacyDefinition(name='type') + ) + conduction_band_min_energy = Quantity( + type=float, + unit="joule", + description=""" + Conduction band minimum energy. + """, + a_legacy=LegacyDefinition(name='conduction_band_min_energy') + ) + valence_band_max_energy = Quantity( + type=float, + unit="joule", + description=""" + Valence band maximum energy. + """, + a_legacy=LegacyDefinition(name='valence_band_max_energy') + ) + conduction_band_min_k_point = Quantity( + type=np.dtype(np.float64), + shape=[3], + unit="1 / meter", + description=""" + Coordinate of the conduction band minimum in k-space. + """, + a_legacy=LegacyDefinition(name='conduction_band_min_k_point') + ) + valence_band_max_k_point = Quantity( + type=np.dtype(np.float64), + shape=[3], + unit="1 / meter", + description=""" + Coordinate of the valence band minimum in k-space. + """, + a_legacy=LegacyDefinition(name='valence_band_max_k_point') + ) + + class section_k_band(MSection): ''' This section stores information on a $k$-band (electronic or vibrational band @@ -2229,6 +2181,61 @@ class section_k_band(MSection): ''', a_legacy=LegacyDefinition(name='band_structure_kind')) + reciprocal_cell = Quantity( + type=np.dtype(np.float64), + shape=[3, 3], + unit="1 / meter", + description=""" + The reciprocal cell within which the band structure is calculated. + """, + a_legacy=LegacyDefinition(name='reciprocal_cell') + ) + + brillouin_zone = Quantity( + type=str, + description=""" + The Brillouin zone that corresponds to the reciprocal cell used in the + band calculation. The Brillouin Zone is defined as a list of vertices + and facets that are encoded with JSON. The vertices are 3D points in + the reciprocal space, and facets are determined by a chain of vertice + indices, with a right-hand ordering determining the surface normal + direction. + { + "vertices": [[3, 2, 1], ...] + "faces": [[0, 1, 2, 3], ...] + } + """, + a_legacy=LegacyDefinition(name='brillouin_zone') + ) + + section_band_gap = SubSection( + sub_section=section_band_gap.m_def, + repeats=False, + a_legacy=LegacyDefinition(name='section_band_gap') + ) + + section_band_gap_spin_up = SubSection( + sub_section=section_band_gap.m_def, + repeats=False, + a_legacy=LegacyDefinition(name='section_band_gap') + ) + + section_band_gap_spin_down = SubSection( + sub_section=section_band_gap.m_def, + repeats=False, + a_legacy=LegacyDefinition(name='section_band_gap') + ) + + is_standard_path = Quantity( + type=bool, + description=""" + Boolean indicating whether the path follows the standard path for this + cell. The AFLOW standard by Setyawan and Curtarolo is used + (https://doi.org/10.1016/j.commatsci.2010.05.010). + """, + a_legacy=LegacyDefinition(name='is_standard_path') + ) + section_k_band_segment = SubSection( sub_section=SectionProxy('section_k_band_segment'), repeats=True, diff --git a/nomad/metainfo/encyclopedia.py b/nomad/metainfo/encyclopedia.py index 1d160896bc950200ad5a02bf8f0aef1451733357..1553793d43bc2ae663d031ea2c66f4d406e8c67c 100644 --- a/nomad/metainfo/encyclopedia.py +++ b/nomad/metainfo/encyclopedia.py @@ -1,6 +1,7 @@ import numpy as np from elasticsearch_dsl import InnerDoc -from nomad.metainfo import MSection, Section, SubSection, Quantity, MEnum, units +from nomad.metainfo import MSection, Section, SubSection, Quantity, MEnum, Reference +from nomad.datamodel.metainfo.public import section_k_band, section_dos class WyckoffVariables(MSection): @@ -408,199 +409,6 @@ class Calculation(MSection): ) -class BandGap(MSection): - m_def = Section( - a_flask=dict(skip_none=True), - a_elastic=dict(type=InnerDoc), - description=""" - Stores information related to a band gap that has bee identified within - the band structure. - """ - ) - value = Quantity( - type=float, - unit=units.J, - description=""" - Band gap energy. - """ - ) - type = Quantity( - type=MEnum("direct", "indirect"), - description=""" - Type of band gap. - """ - ) - conduction_band_min_energy = Quantity( - type=float, - unit=units.J, - description=""" - Conduction band minimum energy. - """ - ) - valence_band_max_energy = Quantity( - type=float, - unit=units.J, - description=""" - Valence band maximum energy. - """ - ) - conduction_band_min_k_point = Quantity( - type=np.dtype(np.float64), - shape=[3], - unit=units.m**(-1), - description=""" - Coordinate of the conduction band minimum in k-space. - """ - ) - valence_band_max_k_point = Quantity( - type=np.dtype(np.float64), - shape=[3], - unit=units.m**(-1), - description=""" - Coordinate of the valence band minimum in k-space. - """ - ) - - -class BandSegment(MSection): - m_def = Section( - a_flask=dict(skip_none=True), - a_elastic=dict(type=InnerDoc), - description=""" - Represents a continuous path segment starting from a specific k-point - and ending in another. - """ - ) - energies = Quantity( - type=np.dtype(np.float64), - shape=["1..2", "1..*", "1..*"], - unit=units.m**(-1), - description=""" - The energies of the bands as a 3D array with size [n_spin_channels, - n_bands, n_k_points]. By default the spin down channel is given first - and the spin up channel is second. - """ - ) - k_points = Quantity( - type=np.dtype(np.float64), - shape=["1..*", 3], - unit=units.m**(-1), - description=""" - Path of the band structure in reciprocal space. - """ - ) - labels = Quantity( - type=np.dtype('str'), - shape=[1], - description=""" - The start end end labels for the k points in this segment. - """ - ) - - -class ElectronicBandStructure(MSection): - m_def = Section( - a_flask=dict(skip_none=True), - a_elastic=dict(type=InnerDoc), - description=""" - Stores information related to an electronic band structure. - """ - ) - scc_index = Quantity( - type=int, - description=""" - Index of the single configuration calculation that contains the band - structure. - """ - ) - fermi_level = Quantity( - type=float, - unit=units.J, - description=""" - Fermi level reported for the band structure. - """ - ) - reciprocal_cell = Quantity( - type=np.dtype(np.float64), - shape=[3, 3], - unit=units.m**(-1), - description=""" - The reciprocal cell within which the band structure is calculated. - """ - ) - brillouin_zone = Quantity( - type=str, - description=""" - The Brillouin zone that corresponds to the reciprocal cell used in the - band calculation. The Brillouin Zone is defined as a list of vertices - and facets that are encoded with JSON. The vertices are 3D points in - the reciprocal space, and facets are determined by a chain of vertice - indices, with a right-hand ordering determining the surface normal - direction. - { - "vertices": [[3, 2, 1], ...] - "faces": [[0, 1, 2, 3], ...] - } - """ - ) - band_gap = SubSection(sub_section=BandGap.m_def, repeats=False) - band_gap_spin_up = SubSection(sub_section=BandGap.m_def, repeats=False) - band_gap_spin_down = SubSection(sub_section=BandGap.m_def, repeats=False) - segments = SubSection(sub_section=BandSegment.m_def, repeats=True) - - is_standard_path = Quantity( - type=bool, - description=""" - Boolean indicating whether the path follows the standard path for this - cell. The AFLOW standard by Setyawan and Curtarolo is used - (https://doi.org/10.1016/j.commatsci.2010.05.010). - """ - ) - - -class ElectronicDOS(MSection): - m_def = Section( - a_flask=dict(skip_none=True), - a_elastic=dict(type=InnerDoc), - description=""" - Stores the electronic density of states (DOS). - """ - ) - scc_index = Quantity( - type=int, - description=""" - Index of the single configuration calculation that contains the density - of states. - """ - ) - fermi_level = Quantity( - type=np.dtype(np.float64), - unit=units.J, - description=""" - Fermi level reported for the density of states. - """ - ) - energies = Quantity( - type=np.dtype(np.float64), - shape=["1..*"], - unit=units.J, - description=""" - Array containing the set of discrete energy values with respect to the - top of the valence band for the density of states (DOS). - """ - ) - values = Quantity( - type=np.dtype(np.float64), - shape=["1..2", "1..*"], - unit=units.J**(-1), - description=""" - Values (number of states for a given energy, the set of discrete energy - values is given in dos_energies) of density of states. The values are - given as states/energy. - """ - ) - - class Properties(MSection): m_def = Section( a_flask=dict(skip_none=True), @@ -612,14 +420,14 @@ class Properties(MSection): ) atomic_density = Quantity( type=float, - unit=units.m**(-3), + unit="1 / m ** 3", description=""" Atomic density of the material (atoms/volume)." """ ) mass_density = Quantity( type=float, - unit=units.kg / units.m**3, + unit="kg / m ** 3", description=""" Mass density of the material. """ @@ -630,8 +438,20 @@ class Properties(MSection): Code dependent energy values, corrected to be per formula unit. """ ) - electronic_band_structure = SubSection(sub_section=ElectronicBandStructure.m_def, repeats=False) - electronic_dos = SubSection(sub_section=ElectronicDOS.m_def, repeats=False) + electronic_band_structure = Quantity( + type=Reference(section_k_band.m_def), + shape=[], + description=""" + Reference to an electronic band structure. + """ + ) + electronic_dos = Quantity( + type=Reference(section_dos.m_def), + shape=[], + description=""" + Reference to an electronic density of states. + """ + ) class section_encyclopedia(MSection): diff --git a/nomad/normalizing/band_structure.py b/nomad/normalizing/band_structure.py index 59a3a82467dbb6baa76c9e438959c7c042c4b272..dec950115d790f57da026b9b65c9f08906ef170f 100644 --- a/nomad/normalizing/band_structure.py +++ b/nomad/normalizing/band_structure.py @@ -119,3 +119,204 @@ try: convert_unit_function("m", "angstrom")(self.cell)) except Exception as e: self.logger.warn("failed to get special points/xtal structure", exc_info=e) + + + # Loop over bands + for band_data in bands: + + # Add reciprocal cell and brillouin zone + self.get_reciprocal_cell(band_data, prim_atoms, orig_atoms) + self.get_brillouin_zone(band_data) + + # If the parser has reported a valence band maximum, we can add + # band gap information. + if valence_band_maximum is not None: + kpoints = [] + energies = [] + segments = band_data.section_k_band_segment + if not segments: + return + + # Loop over segments + for segment_src in segments: + try: + seg_k_points = segment_src.band_k_points + seg_energies = segment_src.band_energies + except Exception: + return + if seg_k_points is None or seg_energies is None: + return + else: + seg_energies = seg_energies.magnitude + + seg_energies = np.swapaxes(seg_energies, 1, 2) + kpoints.append(seg_k_points) + energies.append(seg_energies) + + # Continue to calculate band gaps and other properties. + kpoints = np.concatenate(kpoints, axis=0) + energies = np.concatenate(energies, axis=2) + + # If the parser has reported a valence band maximum, we can add + # band gap information. + self.get_band_gaps(band_data, energies, kpoints, valence_band_maximum) + +def get_band_gaps(self, band_structure: ElectronicBandStructure, energies: np.array, path: np.array) -> None: + """Given the band structure and fermi level, calculates the band gap + for spin channels and also reports the total band gap as the minum gap + found. + """ + # Handle spin channels separately to find gaps for spin up and down + reciprocal_cell = band_structure.reciprocal_cell.magnitude + fermi_level = band_structure.fermi_level + n_channels = energies.shape[0] + + gaps: List[BandGap] = [None, None] + for channel in range(n_channels): + channel_energies = energies[channel, :, :] + lower_defined = False + upper_defined = False + num_bands = channel_energies.shape[0] + band_indices = np.arange(num_bands) + band_minima_idx = channel_energies.argmin(axis=1) + band_maxima_idx = channel_energies.argmax(axis=1) + band_minima = channel_energies[band_indices, band_minima_idx] + band_maxima = channel_energies[band_indices, band_maxima_idx] + + # Add a tolerance to minima and maxima + band_minima_tol = band_minima + config.normalize.fermi_level_precision + band_maxima_tol = band_maxima - config.normalize.fermi_level_precision + + for band_idx in range(num_bands): + band_min = band_minima[band_idx] + band_max = band_maxima[band_idx] + band_min_tol = band_minima_tol[band_idx] + band_max_tol = band_maxima_tol[band_idx] + + # If any of the bands band crosses the Fermi level, there is no + # band gap + if band_min_tol <= fermi_level and band_max_tol >= fermi_level: + break + # Whole band below Fermi level, save the current highest + # occupied band point + elif band_min_tol <= fermi_level and band_max_tol <= fermi_level: + gap_lower_energy = band_max + gap_lower_idx = band_maxima_idx[band_idx] + lower_defined = True + # Whole band above Fermi level, save the current lowest + # unoccupied band point + elif band_min_tol >= fermi_level: + gap_upper_energy = band_min + gap_upper_idx = band_minima_idx[band_idx] + upper_defined = True + break + + # If a highest point of the valence band and a lowest point of the + # conduction band are found, and the difference between them is + # positive, save the information location and value. + if lower_defined and upper_defined and gap_upper_energy - gap_lower_energy >= 0: + + # See if the gap is direct or indirect by comparing the k-point + # locations with some tolerance + k_point_lower = path[gap_lower_idx] + k_point_upper = path[gap_upper_idx] + k_point_distance = self.get_k_space_distance(reciprocal_cell, k_point_lower, k_point_upper) + is_direct_gap = k_point_distance <= config.normalize.k_space_precision + + gap = BandGap() + gap.type = "direct" if is_direct_gap else "indirect" + gap.value = float(gap_upper_energy - gap_lower_energy) + gap.conduction_band_min_k_point = k_point_upper + gap.conduction_band_min_energy = float(gap_upper_energy) + gap.valence_band_max_k_point = k_point_lower + gap.valence_band_max_energy = float(gap_lower_energy) + gaps[channel] = gap + + # For unpolarized calculations we simply report the gap if it is found. + if n_channels == 1: + if gaps[0] is not None: + band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[0]) + # For polarized calculations we report the gap separately for both + # channels. Also we report the smaller gap as the the total gap for the + # calculation. + elif n_channels == 2: + if gaps[0] is not None: + band_structure.m_add_sub_section(ElectronicBandStructure.band_gap_spin_up, gaps[0]) + if gaps[1] is not None: + band_structure.m_add_sub_section(ElectronicBandStructure.band_gap_spin_down, gaps[1]) + if gaps[0] is not None and gaps[1] is not None: + if gaps[0].value <= gaps[1].value: + band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[0]) + else: + band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[1]) + else: + if gaps[0] is not None: + band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[0]) + elif gaps[1] is not None: + band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[1]) + + def get_reciprocal_cell(self, band_structure: ElectronicBandStructure, prim_atoms: Atoms, orig_atoms: Atoms): + """A reciprocal cell for this calculation. If the original unit cell is + not a primitive one, then we will use the one given by spglib. + + If the used code is calculating it's own primitive cell, then the + reciprocal cell used in e.g. band structure calculation might differ + from the one given by spglib. To overcome this we would need to get the + exact reciprocal cell used by the calculation or then test for + different primitive cells + (https://atztogo.github.io/spglib/definition.html#transformation-to-the-primitive-cell) + whether the k-point path stays inside the first Brillouin zone. + + Args: + prim_atoms: The primitive system as detected by MatID. + orig_atoms: The original simulation cell. + + Returns: + Reciprocal cell as 3x3 numpy array. Units are in SI (1/m). + """ + primitive_cell = prim_atoms.get_cell() + source_cell = orig_atoms.get_cell() + + volume_primitive = primitive_cell.volume + volume_source = source_cell.volume + volume_diff = abs(volume_primitive - volume_source) + + if volume_diff > (0.001)**3: + recip_cell = primitive_cell.reciprocal() * 1e10 + else: + recip_cell = source_cell.reciprocal() * 1e10 + + band_structure.reciprocal_cell = recip_cell + + + def get_k_space_distance(self, reciprocal_cell: np.array, point1: np.array, point2: np.array) -> float: + """Used to calculate the Euclidean distance of two points in k-space, + given relative positions in the reciprocal cell. + + Args: + reciprocal_cell: Reciprocal cell. + point1: The first position in k-space. + point2: The second position in k-space. + + Returns: + float: Euclidian distance of the two points in k-space in SI units. + """ + k_point_displacement = np.dot(reciprocal_cell, point1 - point2) + k_point_distance = np.linalg.norm(k_point_displacement) + + return k_point_distance + + def get_brillouin_zone(self, band_structure: ElectronicBandStructure) -> None: + """Returns a dictionary containing the information needed to display + the Brillouin zone for this material. This functionality could be put + into the GUI directly, with the Brillouin zone construction performed + from the reciprocal cell. + + The Brillouin Zone is a Wigner-Seitz cell, and is thus uniquely + defined. It's shape does not depend on the used primitive cell. + """ + recip_cell = band_structure.reciprocal_cell.magnitude + brillouin_zone = atomutils.get_brillouin_zone(recip_cell) + bz_json = json.dumps(brillouin_zone) + band_structure.brillouin_zone = bz_json + diff --git a/nomad/normalizing/encyclopedia/properties.py b/nomad/normalizing/encyclopedia/properties.py index aa3576db266c7c3b9e55f8af1660c0c28b326ff1..d13ebebd120d6e5e9572ec18cfdd051afec0aecb 100644 --- a/nomad/normalizing/encyclopedia/properties.py +++ b/nomad/normalizing/encyclopedia/properties.py @@ -12,25 +12,16 @@ # See the License for the specific language governing permissions and # limitations under the License. -from typing import List import json -import numpy as np -from ase import Atoms from nomad.metainfo.encyclopedia import ( Calculation, Properties, Material, - ElectronicBandStructure, - BandGap, ) from nomad.parsing.legacy import Backend from nomad.metainfo import Section from nomad.normalizing.encyclopedia.context import Context -from nomad import atomutils -from nomad import config - -J_to_Ry = 4.587425e+17 class PropertiesNormalizer(): @@ -41,170 +32,9 @@ class PropertiesNormalizer(): self.backend = backend self.logger = logger - def get_reciprocal_cell(self, band_structure: ElectronicBandStructure, prim_atoms: Atoms, orig_atoms: Atoms): - """A reciprocal cell for this calculation. If the original unit cell is - not a primitive one, then we will use the one given by spglib. - - If the used code is calculating it's own primitive cell, then the - reciprocal cell used in e.g. band structure calculation might differ - from the one given by spglib. To overcome this we would need to get the - exact reciprocal cell used by the calculation or then test for - different primitive cells - (https://atztogo.github.io/spglib/definition.html#transformation-to-the-primitive-cell) - whether the k-point path stays inside the first Brillouin zone. - - Args: - prim_atoms: The primitive system as detected by MatID. - orig_atoms: The original simulation cell. - - Returns: - Reciprocal cell as 3x3 numpy array. Units are in SI (1/m). - """ - primitive_cell = prim_atoms.get_cell() - source_cell = orig_atoms.get_cell() - - volume_primitive = primitive_cell.volume - volume_source = source_cell.volume - volume_diff = abs(volume_primitive - volume_source) - - if volume_diff > (0.001)**3: - recip_cell = primitive_cell.reciprocal() * 1e10 - else: - recip_cell = source_cell.reciprocal() * 1e10 - - band_structure.reciprocal_cell = recip_cell - - def get_band_gaps(self, band_structure: ElectronicBandStructure, energies: np.array, path: np.array) -> None: - """Given the band structure and fermi level, calculates the band gap - for spin channels and also reports the total band gap as the minum gap - found. - """ - # Handle spin channels separately to find gaps for spin up and down - reciprocal_cell = band_structure.reciprocal_cell.magnitude - fermi_level = band_structure.fermi_level - n_channels = energies.shape[0] - - gaps: List[BandGap] = [None, None] - for channel in range(n_channels): - channel_energies = energies[channel, :, :] - lower_defined = False - upper_defined = False - num_bands = channel_energies.shape[0] - band_indices = np.arange(num_bands) - band_minima_idx = channel_energies.argmin(axis=1) - band_maxima_idx = channel_energies.argmax(axis=1) - band_minima = channel_energies[band_indices, band_minima_idx] - band_maxima = channel_energies[band_indices, band_maxima_idx] - - # Add a tolerance to minima and maxima - band_minima_tol = band_minima + config.normalize.fermi_level_precision - band_maxima_tol = band_maxima - config.normalize.fermi_level_precision - - for band_idx in range(num_bands): - band_min = band_minima[band_idx] - band_max = band_maxima[band_idx] - band_min_tol = band_minima_tol[band_idx] - band_max_tol = band_maxima_tol[band_idx] - - # If any of the bands band crosses the Fermi level, there is no - # band gap - if band_min_tol <= fermi_level and band_max_tol >= fermi_level: - break - # Whole band below Fermi level, save the current highest - # occupied band point - elif band_min_tol <= fermi_level and band_max_tol <= fermi_level: - gap_lower_energy = band_max - gap_lower_idx = band_maxima_idx[band_idx] - lower_defined = True - # Whole band above Fermi level, save the current lowest - # unoccupied band point - elif band_min_tol >= fermi_level: - gap_upper_energy = band_min - gap_upper_idx = band_minima_idx[band_idx] - upper_defined = True - break - - # If a highest point of the valence band and a lowest point of the - # conduction band are found, and the difference between them is - # positive, save the information location and value. - if lower_defined and upper_defined and gap_upper_energy - gap_lower_energy >= 0: - - # See if the gap is direct or indirect by comparing the k-point - # locations with some tolerance - k_point_lower = path[gap_lower_idx] - k_point_upper = path[gap_upper_idx] - k_point_distance = self.get_k_space_distance(reciprocal_cell, k_point_lower, k_point_upper) - is_direct_gap = k_point_distance <= config.normalize.k_space_precision - - gap = BandGap() - gap.type = "direct" if is_direct_gap else "indirect" - gap.value = float(gap_upper_energy - gap_lower_energy) - gap.conduction_band_min_k_point = k_point_upper - gap.conduction_band_min_energy = float(gap_upper_energy) - gap.valence_band_max_k_point = k_point_lower - gap.valence_band_max_energy = float(gap_lower_energy) - gaps[channel] = gap - - # For unpolarized calculations we simply report the gap if it is found. - if n_channels == 1: - if gaps[0] is not None: - band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[0]) - # For polarized calculations we report the gap separately for both - # channels. Also we report the smaller gap as the the total gap for the - # calculation. - elif n_channels == 2: - if gaps[0] is not None: - band_structure.m_add_sub_section(ElectronicBandStructure.band_gap_spin_up, gaps[0]) - if gaps[1] is not None: - band_structure.m_add_sub_section(ElectronicBandStructure.band_gap_spin_down, gaps[1]) - if gaps[0] is not None and gaps[1] is not None: - if gaps[0].value <= gaps[1].value: - band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[0]) - else: - band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[1]) - else: - if gaps[0] is not None: - band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[0]) - elif gaps[1] is not None: - band_structure.m_add_sub_section(ElectronicBandStructure.band_gap, gaps[1]) - - def get_k_space_distance(self, reciprocal_cell: np.array, point1: np.array, point2: np.array) -> float: - """Used to calculate the Euclidean distance of two points in k-space, - given relative positions in the reciprocal cell. - - Args: - reciprocal_cell: Reciprocal cell. - point1: The first position in k-space. - point2: The second position in k-space. - - Returns: - float: Euclidian distance of the two points in k-space in SI units. - """ - k_point_displacement = np.dot(reciprocal_cell, point1 - point2) - k_point_distance = np.linalg.norm(k_point_displacement) - - return k_point_distance - - def get_brillouin_zone(self, band_structure: ElectronicBandStructure) -> None: - """Returns a dictionary containing the information needed to display - the Brillouin zone for this material. This functionality could be put - into the GUI directly, with the Brillouin zone construction performed - from the reciprocal cell. - - The Brillouin Zone is a Wigner-Seitz cell, and is thus uniquely - defined. It's shape does not depend on the used primitive cell. - """ - recip_cell = band_structure.reciprocal_cell.magnitude - brillouin_zone = atomutils.get_brillouin_zone(recip_cell) - bz_json = json.dumps(brillouin_zone) - band_structure.brillouin_zone = bz_json - - def band_structure(self, properties: Properties, calc_type: str, material_type: str, context: Context, sec_system: Section) -> None: - """Band structure data following arbitrary path. - - Currently this function is only taking into account the normalized band - structures, and if they have a k-point that is labeled as '?', the - whole band strcture will be ignored. + 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. """ # Band structure data is extracted only from bulk calculations. This is # because for now we need the reciprocal cell of the primitive cell @@ -214,64 +44,42 @@ class PropertiesNormalizer(): if calc_type != Calculation.calculation_type.type.single_point or material_type != Material.material_type.type.bulk: return + # Try to find a band structure. representative_scc = context.representative_scc - orig_atoms = sec_system.m_cache["representative_atoms"] - symmetry_analyzer = sec_system.section_symmetry[0].m_cache["symmetry_analyzer"] - prim_atoms = symmetry_analyzer.get_primitive_system() - - # Try to find an SCC with band structure data, give priority to - # normalized data - for src_name in ["section_k_band_normalized", "section_k_band"]: - bands = getattr(representative_scc, src_name) - if bands is None: - continue - norm = "_normalized" if src_name == "section_k_band_normalized" else "" - - # Loop over bands - for band_data in bands: - band_structure = ElectronicBandStructure() - band_structure.scc_index = int(context.representative_scc_idx) - kpoints = [] - energies = [] - segments = band_data['section_k_band_segment' + norm] - if not segments: - return + bands = representative_scc.section_k_band + if bands is None or len(bands) == 0: + return - # Loop over segments - for segment_src in segments: - try: - seg_k_points = segment_src["band_k_points" + norm] - seg_energies = segment_src["band_energies" + norm] - seg_labels = segment_src['band_segm_labels' + norm] - except Exception: - return - if seg_k_points is None or seg_energies is None or seg_labels is None: - return - else: - seg_energies = seg_energies.magnitude - if "?" in seg_labels: - return + # Loop over bands + representative_band = None + for band in bands: - seg_energies = np.swapaxes(seg_energies, 1, 2) - kpoints.append(seg_k_points) - energies.append(seg_energies) + # Check segments + segments = band.section_k_band_segment + if segments is None or len(segments) == 0: + continue - # Continue to calculate band gaps and other properties. - kpoints = np.concatenate(kpoints, axis=0) - energies = np.concatenate(energies, axis=2) - self.get_reciprocal_cell(band_structure, prim_atoms, orig_atoms) - self.get_brillouin_zone(band_structure) + # 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 - # If we are using a normalized band structure (or the Fermi level - # is defined somehow else), we can add band gap information - fermi_level = None - if src_name == "section_k_band_normalized": - fermi_level = 0.0 - if fermi_level is not None: - band_structure.fermi_level = fermi_level - self.get_band_gaps(band_structure, energies, kpoints) + # 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 - properties.m_add_sub_section(Properties.electronic_band_structure, band_structure) + if representative_band is not None: + properties.electronic_band_structure = representative_band def elastic_constants_matrix(self) -> None: pass @@ -342,5 +150,5 @@ class PropertiesNormalizer(): gcd = context.greatest_common_divisor # Save metainfo - self.band_structure(properties, calc_type, material_type, context, sec_system) + self.electronic_band_structure(properties, calc_type, material_type, context, sec_system) self.energies(properties, gcd, representative_scc)