From 77cfeec98f1b5db36b4bffda29fa506c90ef6d58 Mon Sep 17 00:00:00 2001 From: jrudz Date: Thu, 13 Jun 2024 10:24:28 +0200 Subject: [PATCH] added missing property files by hand --- properties/band_gap.py | 159 +++++++++ properties/band_structure.py | 328 +++++++++++++++++ properties/energies.py | 79 +++++ properties/fermi_surface.py | 53 +++ properties/hopping_matrix.py | 105 ++++++ properties/permittivity.py | 119 +++++++ properties/spectral_profile.py | 619 +++++++++++++++++++++++++++++++++ 7 files changed, 1462 insertions(+) create mode 100644 properties/band_gap.py create mode 100644 properties/band_structure.py create mode 100644 properties/fermi_surface.py create mode 100644 properties/hopping_matrix.py create mode 100644 properties/permittivity.py create mode 100644 properties/spectral_profile.py diff --git a/properties/band_gap.py b/properties/band_gap.py new file mode 100644 index 00000000..c6922a4f --- /dev/null +++ b/properties/band_gap.py @@ -0,0 +1,159 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# 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 +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# 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. +# + +import numpy as np +from structlog.stdlib import BoundLogger +import pint +from typing import Optional + +from nomad.metainfo import Quantity, MEnum, Section, Context + +from nomad_simulations.physical_property import PhysicalProperty + + +class ElectronicBandGap(PhysicalProperty): + """ + Energy difference between the highest occupied electronic state and the lowest unoccupied electronic state. + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicBandGap' + + type = Quantity( + type=MEnum('direct', 'indirect'), + description=""" + Type categorization of the electronic band gap. This quantity is directly related with `momentum_transfer` as by + definition, the electronic band gap is `'direct'` for zero momentum transfer (or if `momentum_transfer` is `None`) and `'indirect'` + for finite momentum transfer. + + Note: in the case of finite `variables`, this quantity refers to all of the `value` in the array. + """, + ) + + momentum_transfer = Quantity( + type=np.float64, + shape=[2, 3], + description=""" + If the electronic band gap is `'indirect'`, the reciprocal momentum transfer for which the band gap is defined + in units of the `reciprocal_lattice_vectors`. The initial and final momentum 3D vectors are given in the first + and second element. Example, the momentum transfer in bulk Si2 happens between the Γ and the (approximately) + X points in the Brillouin zone; thus: + `momentum_transfer = [[0, 0, 0], [0.5, 0.5, 0]]`. + + Note: this quantity only refers to scalar `value`, not to arrays of `value`. + """, + ) + + spin_channel = Quantity( + type=np.int32, + description=""" + Spin channel of the corresponding electronic band gap. It can take values of 0 or 1. + """, + ) + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the electronic band gap. This value has to be positive, otherwise it will + prop an error and be set to None by the `normalize()` function. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.name = self.m_def.name + self.rank = [] + + def validate_values(self, logger: BoundLogger) -> Optional[pint.Quantity]: + """ + Validate the electronic band gap `value` by checking if they are negative and sets them to None if they are. + + Args: + logger (BoundLogger): The logger to log messages. + """ + value = self.value.magnitude + if not isinstance(self.value.magnitude, np.ndarray): # for scalars + value = np.array( + [value] + ) # ! check this when talking with Lauri and Theodore + + # Set the value to 0 when it is negative + if (value < 0).any(): + logger.error('The electronic band gap cannot be defined negative.') + return None + + if not isinstance(self.value.magnitude, np.ndarray): # for scalars + value = value[0] + return value * self.value.u + + def resolve_type(self, logger: BoundLogger) -> Optional[str]: + """ + Resolves the `type` of the electronic band gap based on the stored `momentum_transfer` values. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[str]): The resolved `type` of the electronic band gap. + """ + mtr = self.momentum_transfer if self.momentum_transfer is not None else [] + + # Check if the `momentum_transfer` is [], and return the type and a warning in the log for `indirect` band gaps + if len(mtr) == 0: + if self.type == 'indirect': + logger.warning( + 'The `momentum_transfer` is not stored for an `indirect` band gap.' + ) + return self.type + + # Check if the `momentum_transfer` has at least two elements, and return None if it does not + if len(mtr) == 1: + logger.warning( + 'The `momentum_transfer` should have at least two elements so that the difference can be calculated and the type of electronic band gap can be resolved.' + ) + return None + + # Resolve `type` from the difference between the initial and final momentum transfer + momentum_difference = np.diff(mtr, axis=0) + if (np.isclose(momentum_difference, np.zeros(3))).all(): + return 'direct' + else: + return 'indirect' + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Checks if the `value` is negative and sets it to None if it is. + self.value = self.validate_values(logger) + if self.value is None: + # ? What about deleting the class if `value` is None? + logger.error('The `value` of the electronic band gap is not stored.') + return + + # Resolve the `type` of the electronic band gap from `momentum_transfer`, ONLY for scalar `value` + if isinstance(self.value.magnitude, np.ndarray): + logger.info( + 'We do not support `type` which describe individual elements in an array `value`.' + ) + else: + self.type = self.resolve_type(logger) + diff --git a/properties/band_structure.py b/properties/band_structure.py new file mode 100644 index 00000000..239503be --- /dev/null +++ b/properties/band_structure.py @@ -0,0 +1,328 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# 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 +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# 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. +# + +import numpy as np +from structlog.stdlib import BoundLogger +from typing import Optional, Tuple, Union +import pint + +from nomad.metainfo import Quantity, Section, Context, SubSection + +from nomad_simulations.numerical_settings import KSpace +from nomad_simulations.physical_property import ( + PhysicalProperty, + validate_quantity_wrt_value, +) +from nomad_simulations.properties import ElectronicBandGap, FermiSurface +from nomad_simulations.utils import get_sibling_section + + +class BaseElectronicEigenvalues(PhysicalProperty): + """ + A base section used to define basic quantities for the `ElectronicEigenvalues` and `ElectronicBandStructure` properties. + """ + + iri = '' + + n_bands = Quantity( + type=np.int32, + description=""" + Number of bands / eigenvalues. + """, + ) + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the electronic eigenvalues. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # ! `n_bands` need to be set up during initialization of the class + self.rank = [int(kwargs.get('n_bands'))] + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class ElectronicEigenvalues(BaseElectronicEigenvalues): + """ """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicEigenvalues' + + spin_channel = Quantity( + type=np.int32, + description=""" + Spin channel of the corresponding electronic eigenvalues. It can take values of 0 or 1. + """, + ) + + occupation = Quantity( + type=np.float64, + shape=['*', 'n_bands'], + description=""" + Occupation of the electronic eigenvalues. This is a number depending whether the `spin_channel` has been set or not. + If `spin_channel` is set, then this number is between 0 and 2, where 0 means that the state is unoccupied and 2 means + that the state is fully occupied; if `spin_channel` is not set, then this number is between 0 and 1. The shape of + this quantity is defined as `[K.n_points, K.dimensionality, n_bands]`, where `K` is a `variable` which can + be `KMesh` or `KLinePath`, depending whether the simulation mapped the whole Brillouin zone or just a specific + path. + """, + ) + + highest_occupied = Quantity( + type=np.float64, + unit='joule', + description=""" + Highest occupied electronic eigenvalue. Together with `lowest_unoccupied`, it defines the + electronic band gap. + """, + ) + + lowest_unoccupied = Quantity( + type=np.float64, + unit='joule', + description=""" + Lowest unoccupied electronic eigenvalue. Together with `highest_occupied`, it defines the + electronic band gap. + """, + ) + + # ? Should we add functionalities to handle min/max of the `value` in some specific cases, .e.g, bands around the Fermi level, + # ? core bands separated by gaps, and equivalently, higher-energy valence bands separated by gaps? + + value_contributions = SubSection( + sub_section=BaseElectronicEigenvalues.m_def, + repeats=True, + description=""" + Contributions to the electronic eigenvalues. Example, in the case of a DFT+GW calculation, the GW eigenvalues + are stored under `value`, and each contribution is identified by `label`: + - `'KS'`: Kohn-Sham contribution. This is also stored in the DFT entry under `ElectronicEigenvalues.value`. + - `'KSxc'`: Diagonal matrix elements of the expectation value of the Kohn-Sahm exchange-correlation potential. + - `'SigX'`: Diagonal matrix elements of the exchange self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`. + - `'SigC'`: Diagonal matrix elements of the correlation self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`. + - `'Zk'`: Quasiparticle renormalization factors contribution. This is also stored in the GW entry under `QuasiparticleWeights.value`. + """, + ) + + reciprocal_cell = Quantity( + type=KSpace.reciprocal_lattice_vectors, + description=""" + Reference to the reciprocal lattice vectors stored under `KSpace`. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.name = self.m_def.name + + @validate_quantity_wrt_value(name='occupation') + def order_eigenvalues(self) -> Union[bool, Tuple[pint.Quantity, np.ndarray]]: + """ + Order the eigenvalues based on the `value` and `occupation`. The return `value` and + `occupation` are flattened. + + Returns: + (Union[bool, Tuple[pint.Quantity, np.ndarray]]): The flattened and sorted `value` and `occupation`. If validation + fails, then it returns `False`. + """ + total_shape = np.prod(self.value.shape) + + # Order the indices in the flattened list of `value` + flattened_value = self.value.reshape(total_shape) + flattened_occupation = self.occupation.reshape(total_shape) + sorted_indices = np.argsort(flattened_value, axis=0) + + sorted_value = ( + np.take_along_axis(flattened_value.magnitude, sorted_indices, axis=0) + * flattened_value.u + ) + sorted_occupation = np.take_along_axis( + flattened_occupation, sorted_indices, axis=0 + ) + self.m_cache['sorted_eigenvalues'] = True + return sorted_value, sorted_occupation + + def resolve_homo_lumo_eigenvalues( + self, + ) -> Tuple[Optional[pint.Quantity], Optional[pint.Quantity]]: + """ + Resolve the `highest_occupied` and `lowest_unoccupied` eigenvalues by performing a binary search on the + flattened and sorted `value` and `occupation`. If these quantities already exist, overwrite them or return + them if it is not possible to resolve from `value` and `occupation`. + + Returns: + (Tuple[Optional[pint.Quantity], Optional[pint.Quantity]]): The `highest_occupied` and + `lowest_unoccupied` eigenvalues. + """ + # Sorting `value` and `occupation` + if not self.order_eigenvalues(): # validation fails + if self.highest_occupied is not None and self.lowest_unoccupied is not None: + return self.highest_occupied, self.lowest_unoccupied + return None, None + sorted_value, sorted_occupation = self.order_eigenvalues() + sorted_value_unit = sorted_value.u + sorted_value = sorted_value.magnitude + + # Binary search ot find the transition point between `occupation = 2` and `occupation = 0` + tolerance = 1e-3 # TODO add tolerance from config fields + homo = self.highest_occupied + lumo = self.lowest_unoccupied + mid = np.searchsorted(sorted_occupation <= tolerance, True) - 1 + if mid >= 0 and mid < len(sorted_occupation) - 1: + if sorted_occupation[mid] > 0 and ( + sorted_occupation[mid + 1] >= -tolerance + and sorted_occupation[mid + 1] <= tolerance + ): + homo = sorted_value[mid] * sorted_value_unit + lumo = sorted_value[mid + 1] * sorted_value_unit + + return homo, lumo + + def extract_band_gap(self) -> Optional[ElectronicBandGap]: + """ + Extract the electronic band gap from the `highest_occupied` and `lowest_unoccupied` eigenvalues. + If the difference of `highest_occupied` and `lowest_unoccupied` is negative, the band gap `value` is set to 0.0. + + Returns: + (Optional[ElectronicBandGap]): The extracted electronic band gap section to be stored in `Outputs`. + """ + band_gap = None + homo, lumo = self.resolve_homo_lumo_eigenvalues() + if homo and lumo: + band_gap = ElectronicBandGap(is_derived=True, physical_property_ref=self) + + if (lumo - homo).magnitude < 0: + band_gap.value = 0.0 + else: + band_gap.value = lumo - homo + return band_gap + + # TODO fix this method once `FermiSurface` property is implemented + def extract_fermi_surface(self, logger: BoundLogger) -> Optional[FermiSurface]: + """ + Extract the Fermi surface for metal systems and using the `FermiLevel.value`. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[FermiSurface]): The extracted Fermi surface section to be stored in `Outputs`. + """ + # Check if the system has a finite band gap + homo, lumo = self.resolve_homo_lumo_eigenvalues() + if (homo and lumo) and (lumo - homo).magnitude > 0: + return None + + # Get the `fermi_level.value` + fermi_level = get_sibling_section( + section=self, sibling_section_name='fermi_level', logger=logger + ) + if fermi_level is None: + logger.warning( + 'Could not extract the `FermiSurface`, because `FermiLevel` is not stored.' + ) + return None + fermi_level_value = fermi_level.value.magnitude + + # Extract values close to the `fermi_level.value` + tolerance = 1e-8 # TODO change this for a config field + fermi_indices = np.logical_and( + self.value.magnitude >= (fermi_level_value - tolerance), + self.value.magnitude <= (fermi_level_value + tolerance), + ) + fermi_values = self.value[fermi_indices] + + # Store `FermiSurface` values + # ! This is wrong (!) the `value` should be the `KMesh.points`, not the `ElectronicEigenvalues.value` + fermi_surface = FermiSurface( + n_bands=self.n_bands, + is_derived=True, + physical_property_ref=self, + ) + fermi_surface.variables = self.variables + fermi_surface.value = fermi_values + return fermi_surface + + def resolve_reciprocal_cell(self) -> Optional[pint.Quantity]: + """ + Resolve the reciprocal cell from the `KSpace` numerical settings section. + + Returns: + Optional[pint.Quantity]: _description_ + """ + numerical_settings = self.m_xpath( + 'm_parent.m_parent.model_method[-1].numerical_settings', dict=False + ) + if numerical_settings is None: + return None + k_space = None + for setting in numerical_settings: + if isinstance(setting, KSpace): + k_space = setting + break + if k_space is None: + return None + return k_space + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Resolve `highest_occupied` and `lowest_unoccupied` eigenvalues + self.highest_occupied, self.lowest_unoccupied = ( + self.resolve_homo_lumo_eigenvalues() + ) + + # `ElectronicBandGap` extraction + band_gap = self.extract_band_gap() + if band_gap is not None: + self.m_parent.electronic_band_gaps.append(band_gap) + + # TODO uncomment once `FermiSurface` property is implemented + # `FermiSurface` extraction + # fermi_surface = self.extract_fermi_surface(logger) + # if fermi_surface is not None: + # self.m_parent.fermi_surfaces.append(fermi_surface) + + # Resolve `reciprocal_cell` from the `KSpace` numerical settings section + self.reciprocal_cell = self.resolve_reciprocal_cell() + + +class ElectronicBandStructure(ElectronicEigenvalues): + """ + Accessible energies by the charges (electrons and holes) in the reciprocal space. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicBandStructure' + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + diff --git a/properties/energies.py b/properties/energies.py index 800f69d4..8ef7c7fb 100644 --- a/properties/energies.py +++ b/properties/energies.py @@ -38,6 +38,85 @@ from nomad.metainfo import Quantity, SubSection, MEnum from .physical_property import PhysicalProperty +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# 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 +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# 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. +# + +import numpy as np + +from nomad.metainfo import Quantity, Section, Context + +from nomad_simulations.physical_property import PhysicalProperty + + +class FermiLevel(PhysicalProperty): + """ + Energy required to add or extract a charge from a material at zero temperature. It can be also defined as the chemical potential at zero temperature. + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + iri = 'http://fairmat-nfdi.eu/taxonomy/FermiLevel' + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the Fermi level. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.rank = [] + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class ChemicalPotential(PhysicalProperty): + """ + Free energy cost of adding or extracting a particle from a thermodynamic system. + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + iri = 'http://fairmat-nfdi.eu/taxonomy/ChemicalPotential' + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the chemical potential. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.rank = [] + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) class EnergyContributions(PhysicalProperty): diff --git a/properties/fermi_surface.py b/properties/fermi_surface.py new file mode 100644 index 00000000..4b473a9a --- /dev/null +++ b/properties/fermi_surface.py @@ -0,0 +1,53 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# 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 +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# 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. +# + +import numpy as np + +from nomad.metainfo import Quantity, Section, Context + +from nomad_simulations.physical_property import PhysicalProperty + + +# TODO This class is not implemented yet. @JosePizarro3 will work in another PR to implement it. +class FermiSurface(PhysicalProperty): + """ + Energy boundary in reciprocal space that separates the filled and empty electronic states in a metal. + It is related with the crossing points in reciprocal space by the chemical potential or, equivalently at + zero temperature, the Fermi level. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/FermiSurface' + + n_bands = Quantity( + type=np.int32, + description=""" + Number of bands / eigenvalues. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # ! `n_bands` need to be set up during initialization of the class + self.rank = [int(kwargs.get('n_bands'))] + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + diff --git a/properties/hopping_matrix.py b/properties/hopping_matrix.py new file mode 100644 index 00000000..cebeea71 --- /dev/null +++ b/properties/hopping_matrix.py @@ -0,0 +1,105 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# 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 +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# 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. +# + +import numpy as np + +from nomad.metainfo import Quantity, Section, Context + +from nomad_simulations.physical_property import PhysicalProperty + + +class HoppingMatrix(PhysicalProperty): + """ + Transition probability between two atomic orbitals in a tight-binding model. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/HoppingMatrix' + + n_orbitals = Quantity( + type=np.int32, + description=""" + Number of orbitals in the tight-binding model. + """, + ) + + degeneracy_factors = Quantity( + type=np.int32, + shape=['*'], + description=""" + Degeneracy of each Wigner-Seitz point. + """, + ) + + value = Quantity( + type=np.complex128, + unit='joule', + description=""" + Value of the hopping matrix in joules. The elements are complex numbers defined for each Wigner-Seitz point and + each pair of orbitals; thus, `rank = [n_orbitals, n_orbitals]`. Note this contains also the onsite values, i.e., + it includes the Wigner-Seitz point (0, 0, 0), hence the `CrystalFieldSplitting` values. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # ! n_orbitals need to be set up during initialization of the class + self.rank = [int(kwargs.get('n_orbitals')), int(kwargs.get('n_orbitals'))] + self.name = self.m_def.name + + # TODO add normalization to extract DOS, band structure, etc, properties from `HoppingMatrix` + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class CrystalFieldSplitting(PhysicalProperty): + """ + Energy difference between the degenerated orbitals of an ion in a crystal field environment. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/CrystalFieldSplitting' + + n_orbitals = Quantity( + type=np.int32, + description=""" + Number of orbitals in the tight-binding model. + """, + ) + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the crystal field splittings in joules. This is the intra-orbital local contribution, i.e., the same orbital + at the same Wigner-Seitz point (0, 0, 0). + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # ! `n_orbitals` need to be set up during initialization of the class + self.rank = [int(kwargs.get('n_orbitals'))] + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + diff --git a/properties/permittivity.py b/properties/permittivity.py new file mode 100644 index 00000000..d97c5cd4 --- /dev/null +++ b/properties/permittivity.py @@ -0,0 +1,119 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# 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 +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# 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. +# + +import numpy as np +from structlog.stdlib import BoundLogger +from typing import Optional, List + +from nomad.metainfo import Quantity, Section, Context, MEnum + +from nomad_simulations.physical_property import PhysicalProperty +from nomad_simulations.utils import get_variables +from nomad_simulations.variables import Frequency, KMesh +from nomad_simulations.properties.spectral_profile import AbsorptionSpectrum + + +# TODO add `DielectricStrength` when we have examples and understand how to extract it from the `Permittivity` tensor. + + +class Permittivity(PhysicalProperty): + """ + Response of the material to polarize in the presence of an electric field. + + Alternative names: `DielectricFunction`. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/Permittivity' + + type = Quantity( + type=MEnum('static', 'dynamic'), + description=""" + Type of permittivity which allows to identify if the permittivity depends on the frequency or not. + """, + ) + + value = Quantity( + type=np.complex128, + # unit='joule', # TODO check units (they have to match `SpectralProfile.value`) + description=""" + Value of the permittivity tensor. If the value does not depend on the scattering vector `q`, then we + can extract the optical absorption spectrum from the imaginary part of the permittivity tensor (this is also called + macroscopic dielectric function). + """, + ) + + # ? We need use cases to understand if we need to define contributions to the permittivity tensor. + # ? `ionic` and `electronic` contributions are common in the literature. + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.rank = [3, 3] + self.name = self.m_def.name + self._axes_map = ['xx', 'yy', 'zz'] + + def resolve_type(self) -> str: + frequencies = get_variables(self.variables, Frequency) + if len(frequencies) > 0: + return 'dynamic' + return 'static' + + def extract_absorption_spectra( + self, logger: BoundLogger + ) -> Optional[List[AbsorptionSpectrum]]: + """ + Extract the absorption spectrum from the imaginary part of the permittivity tensor. + """ + # If the `pemittivity` depends on the scattering vector `q`, then we cannot extract the absorption spectrum + q_mesh = get_variables(self.variables, KMesh) + if len(q_mesh) > 0: + logger.warning( + 'The `permittivity` depends on the scattering vector `q`, so that we cannot extract the absorption spectrum.' + ) + return None + # Extract the `Frequency` variable to extract the absorption spectrum + frequencies = get_variables(self.variables, Frequency) + if len(frequencies) == 0: + logger.warning( + 'The `permittivity` does not have a `Frequency` variable to extract the absorption spectrum.' + ) + return None + # Define the `absorption_spectra` for each principal direction along the diagonal of the `Permittivity.value` as the imaginary part + spectra = [] + for i in range(3): + val = self.value[:, i, i].imag + absorption_spectrum = AbsorptionSpectrum( + axis=self._axes_map[i], variables=frequencies + ) + absorption_spectrum.value = val + absorption_spectrum.physical_property_ref = self + spectra.append(absorption_spectrum) + return spectra + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Resolve the `type` of permittivity + self.type = self.resolve_type() + + # `AbsorptionSpectrum` extraction + absorption_spectra = self.extract_absorption_spectra(logger) + if absorption_spectra is not None: + self.m_parent.absorption_spectrum = absorption_spectra + diff --git a/properties/spectral_profile.py b/properties/spectral_profile.py new file mode 100644 index 00000000..f2c894f0 --- /dev/null +++ b/properties/spectral_profile.py @@ -0,0 +1,619 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# 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 +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# 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. +# + +import numpy as np +from structlog.stdlib import BoundLogger +from typing import Optional, List, Dict +import pint + +from nomad import config +from nomad.metainfo import Quantity, SubSection, Section, Context, MEnum + +from nomad_simulations.utils import get_sibling_section, get_variables +from nomad_simulations.physical_property import PhysicalProperty +from nomad_simulations.variables import Energy2 as Energy +from nomad_simulations.atoms_state import AtomsState, OrbitalsState +from nomad_simulations.properties.band_gap import ElectronicBandGap + + +class SpectralProfile(PhysicalProperty): + """ + A base section used to define the spectral profile. + """ + + value = Quantity( + type=np.float64, + description=""" + The value of the intensities of a spectral profile in arbitrary units. + """, + ) # TODO check units and normalization_factor of DOS and Spectras and see whether they can be merged + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.rank = [] + + def is_valid_spectral_profile(self) -> bool: + """ + Check if the spectral profile is valid, i.e., if all `value` are defined positive. + + Returns: + (bool): True if the spectral profile is valid, False otherwise. + """ + if (self.value < 0.0).any(): + return False + return True + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + if self.is_valid_spectral_profile() is False: + logger.error( + 'Invalid negative intensities found: could not validate spectral profile.' + ) + return + + +class DOSProfile(SpectralProfile): + """ + A base section used to define the `value` of the `ElectronicDensityOfState` property. This is useful when containing + contributions for `projected_dos` with the correct unit. + """ + + value = Quantity( + type=np.float64, + unit='1/joule', + description=""" + The value of the electronic DOS. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + + def resolve_pdos_name(self, logger: BoundLogger) -> Optional[str]: + """ + Resolve the `name` of the projected `DOSProfile` from the `entity_ref` section. This is resolved as: + - `'atom X'` with 'X' being the chemical symbol for `AtomsState` references. + - `'orbital Y X'` with 'X' being the chemical symbol and 'Y' the orbital label for `OrbitalsState` references. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[str]): The resolved `name` of the projected DOS profile. + """ + if self.entity_ref is None: + logger.warning( + 'The `entity_ref` is not set for the DOS profile. Could not resolve the `name`.' + ) + return None + + # Resolve the `name` from the `entity_ref` + name = None + if isinstance(self.entity_ref, AtomsState): + name = f'atom {self.entity_ref.chemical_symbol}' + elif isinstance(self.entity_ref, OrbitalsState): + name = f'orbital {self.entity_ref.l_quantum_symbol}{self.entity_ref.ml_quantum_symbol} {self.entity_ref.m_parent.chemical_symbol}' + return name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # We resolve + self.name = self.resolve_pdos_name(logger) + + +class ElectronicDensityOfStates(DOSProfile): + """ + Number of electronic states accessible for the charges per energy and per volume. + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicDensityOfStates' + + spin_channel = Quantity( + type=np.int32, + description=""" + Spin channel of the corresponding electronic DOS. It can take values of 0 or 1. + """, + ) + + # TODO clarify the role of `energies_origin` once `ElectronicEigenvalues` is implemented + energies_origin = Quantity( + type=np.float64, + unit='joule', + description=""" + Energy level denoting the origin along the energy axis, used for comparison and visualization. It is + defined as the `ElectronicEigenvalues.highest_occupied_energy`. + """, + ) + + normalization_factor = Quantity( + type=np.float64, + description=""" + Normalization factor for electronic DOS to get a cell-independent intensive DOS. The cell-independent + intensive DOS is as the integral from the lowest (most negative) energy to the Fermi level for a neutrally + charged system (i.e., the sum of `AtomsState.charge` is zero). + """, + ) + + # ? Do we want to store the integrated value here os as part of an nomad-analysis tool? Check `dos_integrator.py` module in dos normalizer repository + # value_integrated = Quantity( + # type=np.float64, + # description=""" + # The cumulative intensities integrated from from the lowest (most negative) energy to the Fermi level. + # """, + # ) + + projected_dos = SubSection( + sub_section=DOSProfile.m_def, + repeats=True, + description=""" + Projected DOS. It can be atom- (different elements in the unit cell) or orbital-projected. These can be calculated in a cascade as: + - If the total DOS is not present, we sum all atom-projected DOS to obtain it. + - If the atom-projected DOS is not present, we sum all orbital-projected DOS to obtain it. + Note: the cover given by summing up contributions is not perfect, and will depend on the projection functions used. + + In `projected_dos`, `name` and `entity_ref` must be set in order for normalization to work: + - The `entity_ref` is the `OrbitalsState` or `AtomsState` sections. + - The `name` of the projected DOS should be `'atom X'` or `'orbital Y X'`, with 'X' being the chemical symbol and 'Y' the orbital label. + These can be extracted from `entity_ref`. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.name = self.m_def.name + + def resolve_energies_origin( + self, + energies: pint.Quantity, + fermi_level: Optional[pint.Quantity], + logger: BoundLogger, + ) -> Optional[pint.Quantity]: + """ + Resolve the origin of reference for the energies from the sibling `ElectronicEigenvalues` section and its + `highest_occupied` level, or if this does not exist, from the `fermi_level` value as extracted from the sibling property, `FermiLevel`. + + Args: + fermi_level (Optional[pint.Quantity]): The resolved Fermi level. + energies (pint.Quantity): The grid points of the `Energy` variable. + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[pint.Quantity]): The resolved origin of reference for the energies. + """ + # Check if the variables contain more than one variable (different than Energy) + # ? Is this correct or should be use the index of energies to extract the proper shape element in `self.value` being used for `dos_values`? + if len(self.variables) > 1: + logger.warning( + 'The ElectronicDensityOfStates section contains more than one variable. We cannot extract the energy reference.' + ) + return None + + # Extract the `ElectronicEigenvalues` section to get the `highest_occupied` and `lowest_unoccupied` energies + # TODO implement once `ElectronicEigenvalues` is in the schema + eigenvalues = get_sibling_section( + section=self, sibling_section_name='electronic_eigenvalues', logger=logger + ) # we consider `index_sibling` to be 0 + highest_occupied_energy = ( + eigenvalues.highest_occupied if eigenvalues is not None else None + ) + lowest_unoccupied_energy = ( + eigenvalues.lowest_unoccupied if eigenvalues is not None else None + ) + # and set defaults for `highest_occupied_energy` and `lowest_unoccupied_energy` in `m_cache` + if highest_occupied_energy is not None: + self.m_cache['highest_occupied_energy'] = highest_occupied_energy + if lowest_unoccupied_energy is not None: + self.m_cache['lowest_unoccupied_energy'] = lowest_unoccupied_energy + + # Set thresholds for the energies and values + energy_threshold = config.normalize.band_structure_energy_tolerance + value_threshold = 1e-8 # The DOS value that is considered to be zero + + # Check that the closest `energies` to the energy reference is not too far away. + # If it is very far away, normalization may be very inaccurate and we do not report it. + dos_values = self.value.magnitude + eref = highest_occupied_energy if fermi_level is None else fermi_level + fermi_idx = (np.abs(energies - eref)).argmin() + fermi_energy_closest = energies[fermi_idx] + distance = np.abs(fermi_energy_closest - eref) + single_peak_fermi = False + if distance.magnitude <= energy_threshold: + # See if there are zero values close below the energy reference. + idx = fermi_idx + idx_descend = fermi_idx + while True: + try: + value = dos_values[idx] + energy_distance = np.abs(eref - energies[idx]) + except IndexError: + break + if energy_distance.magnitude > energy_threshold: + break + if value <= value_threshold: + idx_descend = idx + break + idx -= 1 + + # See if there are zero values close above the fermi energy. + idx = fermi_idx + idx_ascend = fermi_idx + while True: + try: + value = dos_values[idx] + energy_distance = np.abs(eref - energies[idx]) + except IndexError: + break + if energy_distance.magnitude > energy_threshold: + break + if value <= value_threshold: + idx_ascend = idx + break + idx += 1 + + # If there is a single peak at fermi energy, no + # search needs to be performed. + if idx_ascend != fermi_idx and idx_descend != fermi_idx: + self.m_cache['highest_occupied_energy'] = fermi_energy_closest + self.m_cache['lowest_unoccupied_energy'] = fermi_energy_closest + single_peak_fermi = True + + if not single_peak_fermi: + # Look for highest occupied energy below the descend index + idx = idx_descend + while True: + try: + value = dos_values[idx] + except IndexError: + break + if value > value_threshold: + idx = idx if idx == idx_descend else idx + 1 + self.m_cache['highest_occupied_energy'] = energies[idx] + break + idx -= 1 + # Look for lowest unoccupied energy above idx_ascend + idx = idx_ascend + while True: + try: + value = dos_values[idx] + except IndexError: + break + if value > value_threshold: + idx = idx if idx == idx_ascend else idx - 1 + self.m_cache['highest_occupied_energy'] = energies[idx] + break + idx += 1 + + # Return the `highest_occupied_energy` as the `energies_origin`, or the `fermi_level` if it is not None + energies_origin = self.m_cache.get('highest_occupied_energy') + if energies_origin is None: + energies_origin = fermi_level + return energies_origin + + def resolve_normalization_factor(self, logger: BoundLogger) -> Optional[float]: + """ + Resolve the `normalization_factor` for the electronic DOS to get a cell-independent intensive DOS. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[float]): The normalization factor. + """ + # Get the `ModelSystem` as referenced in the `Outputs.model_system_ref` + model_system = get_sibling_section( + section=self, sibling_section_name='model_system_ref', logger=logger + ) + if model_system is None: + logger.warning( + 'Could not resolve the referenced `ModelSystem` in the `Outputs`.' + ) + return None + + # Get the originally parsed `AtomicCell`, which is the first element stored in `ModelSystem.cell` of name `'AtomicCell'` + atomic_cell = None + for cell in model_system.cell: + if cell.name == 'AtomicCell': # we get the originally parsed `AtomicCell` + atomic_cell = cell + break + if atomic_cell is None: + logger.warning( + 'Could not resolve the `AtomicCell` from the referenced `ModelSystem`.' + ) + return None + + # Get the `atoms_state` and their `atomic_number` from the `AtomicCell` + if atomic_cell.atoms_state is None or len(atomic_cell.atoms_state) == 0: + logger.warning('Could not resolve the `atoms_state` from the `AtomicCell`.') + return None + atomic_numbers = [atom.atomic_number for atom in atomic_cell.atoms_state] + + # Return `normalization_factor` depending if the calculation is spin polarized or not + if self.spin_channel is not None: + normalization_factor = 1 / (2 * sum(atomic_numbers)) + else: + normalization_factor = 1 / sum(atomic_numbers) + return normalization_factor + + def extract_band_gap(self) -> Optional[ElectronicBandGap]: + """ + Extract the electronic band gap from the `highest_occupied_energy` and `lowest_unoccupied_energy` stored + in `m_cache` from `resolve_energies_origin()`. If the difference of `highest_occupied_energy` and + `lowest_unoccupied_energy` is negative, the band gap `value` is set to 0.0. + + Returns: + (Optional[ElectronicBandGap]): The extracted electronic band gap section to be stored in `Outputs`. + """ + band_gap = None + homo = self.m_cache.get('highest_occupied_energy') + lumo = self.m_cache.get('lowest_unoccupied_energy') + if homo and lumo: + band_gap = ElectronicBandGap() + band_gap.is_derived = True + band_gap.physical_property_ref = self + + if (homo - lumo).magnitude < 0: + band_gap.value = 0.0 + else: + band_gap.value = homo - lumo + return band_gap + + def extract_projected_dos( + self, type: str, logger: BoundLogger + ) -> List[Optional[DOSProfile]]: + """ + Extract the projected DOS from the `projected_dos` section and the specified `type`. + + Args: + type (str): The type of the projected DOS to extract. It can be `'atom'` or `'orbital'`. + + Returns: + (DOSProfile): The extracted projected DOS. + """ + extracted_pdos = [] + for pdos in self.projected_dos: + # We make sure each PDOS is normalized + pdos.normalize(None, logger) + + # Initial check for `name` and `entity_ref` + if ( + pdos.name is None + or pdos.entity_ref is None + or len(pdos.entity_ref) == 0 + ): + logger.warning( + '`name` or `entity_ref` are not set for `projected_dos` and they are required for normalization to work.' + ) + return None + + if type in pdos.name: + extracted_pdos.append(pdos) + return extracted_pdos + + def generate_from_projected_dos( + self, logger: BoundLogger + ) -> Optional[pint.Quantity]: + """ + Generate the total `value` of the electronic DOS from the `projected_dos` contributions. If the `projected_dos` + is not present, it returns `None`. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[pint.Quantity]): The total `value` of the electronic DOS. + """ + if self.projected_dos is None or len(self.projected_dos) == 0: + return None + + # Extract `Energy` variables + energies = get_variables(self.variables, Energy) + if len(energies) != 1: + logger.warning( + 'The `ElectronicDensityOfStates` does not contain an `Energy` variable to extract the DOS.' + ) + return None + + # We distinguish between orbital and atom `projected_dos` + orbital_projected = self.extract_projected_dos('orbital', logger) + atom_projected = self.extract_projected_dos('atom', logger) + + # Extract `atom_projected` from `orbital_projected` by summing up the `orbital_projected` contributions for each atom + if len(atom_projected) == 0: + atom_data: Dict[AtomsState, List[DOSProfile]] = {} + for orb_pdos in orbital_projected: + # `entity_ref` is the `OrbitalsState` section, whose parent is `AtomsState` + entity_ref = orb_pdos.entity_ref.m_parent + if entity_ref in atom_data: + atom_data[entity_ref].append(orb_pdos) + else: + atom_data[entity_ref] = [orb_pdos] + for ref, data in atom_data.items(): + atom_dos = DOSProfile( + name=f'atom {ref.chemical_symbol}', + entity_ref=ref, + variables=energies, + ) + orbital_values = [ + dos.value.magnitude for dos in data + ] # to avoid warnings from pint + orbital_unit = data[0].value.u + atom_dos.value = np.sum(orbital_values, axis=0) * orbital_unit + atom_projected.append(atom_dos) + # We concatenate the `atom_projected` to the `projected_dos` + self.projected_dos = orbital_projected + atom_projected + + # Extract `value` from `atom_projected` by summing up the `atom_projected` contributions + value = self.value + if value is None: + atom_values = [ + dos.value.magnitude for dos in atom_projected + ] # to avoid warnings from pint + atom_unit = atom_projected[0].value.u + value = np.sum(atom_values, axis=0) * atom_unit + return value + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Initial check to see if `variables` contains the required `Energy` variable + energies = get_variables(self.variables, Energy) + if len(energies) != 0: + return + energies = energies[0].points + + # Resolve `fermi_level` from a sibling section with respect to `ElectronicDensityOfStates` + fermi_level = get_sibling_section( + section=self, sibling_section_name='fermi_level', logger=logger + ) # * we consider `index_sibling` to be 0 + if fermi_level is not None: + fermi_level = fermi_level.value + # and the `energies_origin` from the sibling `ElectronicEigenvalues` section + self.energies_origin = self.resolve_energies_origin( + energies, fermi_level, logger + ) + if self.energies_origin is None: + logger.info('Could not resolve the `energies_origin` for the DOS') + + # Resolve `normalization_factor` + if self.normalization_factor is None: + self.normalization_factor = self.resolve_normalization_factor(logger) + + # `ElectronicBandGap` extraction + band_gap = self.extract_band_gap() + if band_gap is not None: + self.m_parent.electronic_band_gap.append(band_gap) + + # Total `value` extraction from `projected_dos` + value_from_pdos = self.generate_from_projected_dos(logger) + if self.value is None and value_from_pdos is not None: + logger.info( + 'The `ElectronicDensityOfStates.value` is not stored. We will attempt to obtain it by summing up projected DOS contributions, if these are present.' + ) + self.value = value_from_pdos + + +class AbsorptionSpectrum(SpectralProfile): + """ """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + axis = Quantity( + type=MEnum('xx', 'yy', 'zz'), + description=""" + Axis of the absorption spectrum. This is related with the polarization direction, and can be seen as the + principal term in the tensor `Permittivity.value` (see permittivity.py module). + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # Set the name of the section + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class XASSpectrum(AbsorptionSpectrum): + """ + X-ray Absorption Spectrum (XAS). + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + xanes_spectrum = SubSection( + sub_section=AbsorptionSpectrum.m_def, + description=""" + X-ray Absorption Near Edge Structure (XANES) spectrum. + """, + repeats=False, + ) + + exafs_spectrum = SubSection( + sub_section=AbsorptionSpectrum.m_def, + description=""" + Extended X-ray Absorption Fine Structure (EXAFS) spectrum. + """, + repeats=False, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # Set the name of the section + self.name = self.m_def.name + + def generate_from_contributions(self, logger: BoundLogger) -> None: + """ + Generate the `value` of the XAS spectrum by concatenating the XANES and EXAFS contributions. It also concatenates + the `Energy` grid points of the XANES and EXAFS parts. + + Args: + logger (BoundLogger): The logger to log messages. + """ + # TODO check if this method is general enough + if self.xanes_spectrum is not None and self.exafs_spectrum is not None: + # Concatenate XANE and EXAFS `Energy` grid points + xanes_variables = get_variables(self.xanes_spectrum.variables, Energy) + exafs_variables = get_variables(self.exafs_spectrum.variables, Energy) + if len(xanes_variables) == 0 or len(exafs_variables) == 0: + logger.warning( + 'Could not extract the `Energy` grid points from XANES or EXAFS.' + ) + return + xanes_energies = xanes_variables[0].points + exafs_energies = exafs_variables[0].points + if xanes_energies.max() > exafs_energies.min(): + logger.warning( + 'The XANES `Energy` grid points are not below the EXAFS `Energy` grid points.' + ) + return + self.variables = [ + Energy(points=np.concatenate([xanes_energies, exafs_energies])) + ] + # Concatenate XANES and EXAFS `value` if they have the same shape ['n_energies'] + try: + self.value = np.concatenate( + [self.xanes_spectrum.value, self.exafs_spectrum.value] + ) + except ValueError: + logger.warning( + 'The XANES and EXAFS `value` have different shapes. Could not concatenate the values.' + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + if self.value is None: + logger.info( + 'The `XASSpectrum.value` is not stored. We will attempt to obtain it by combining the XANES and EXAFS parts if these are present.' + ) + self.generate_from_contributions(logger) +