diff --git a/.gitignore b/.gitignore index e818705e..dd3f658c 100644 --- a/.gitignore +++ b/.gitignore @@ -122,6 +122,7 @@ celerybeat.pid # Environments .env .venv +.pyenv env/ venv/ ENV/ diff --git a/src/nomad_simulations/atoms_state.py b/src/nomad_simulations/atoms_state.py index b15f3b8d..c9ff423f 100644 --- a/src/nomad_simulations/atoms_state.py +++ b/src/nomad_simulations/atoms_state.py @@ -25,12 +25,13 @@ from nomad.metainfo import Quantity, SubSection, MEnum, Section, Context from nomad.datamodel.data import ArchiveSection +from nomad.datamodel.metainfo.basesections import Entity from nomad.datamodel.metainfo.annotations import ELNAnnotation from .utils import RussellSaundersState -class OrbitalsState(ArchiveSection): +class OrbitalsState(Entity): """ A base section used to define the orbital state of an atom. """ @@ -537,7 +538,7 @@ def normalize(self, archive, logger) -> None: ) -class AtomsState(ArchiveSection): +class AtomsState(Entity): """ A base section to define each atom state information. """ diff --git a/src/nomad_simulations/outputs.py b/src/nomad_simulations/outputs.py index ff154a66..7d4bfec1 100644 --- a/src/nomad_simulations/outputs.py +++ b/src/nomad_simulations/outputs.py @@ -28,6 +28,10 @@ from .numerical_settings import SelfConsistency from .properties import ( ElectronicBandGap, + FermiLevel, + ChemicalPotential, + ElectronicDensityOfStates, + XASSpectra, ) @@ -53,20 +57,22 @@ class Outputs(ArchiveSection): a_eln=ELNAnnotation(component='ReferenceEditQuantity'), ) - custom_physical_property = SubSection( - sub_section=PhysicalProperty.m_def, - repeats=True, - description=""" - A custom physical property used to store properties not yet covered by the NOMAD schema. - """, - ) - # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # List of properties # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # + fermi_level = SubSection(sub_section=FermiLevel.m_def, repeats=True) + + chemical_potential = SubSection(sub_section=ChemicalPotential.m_def, repeats=True) + electronic_band_gap = SubSection(sub_section=ElectronicBandGap.m_def, repeats=True) + electronic_dos = SubSection( + sub_section=ElectronicDensityOfStates.m_def, repeats=True + ) + + xas_spectra = SubSection(sub_section=XASSpectra.m_def, repeats=True) + # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # diff --git a/src/nomad_simulations/physical_property.py b/src/nomad_simulations/physical_property.py index b013fa5b..afb0fab8 100644 --- a/src/nomad_simulations/physical_property.py +++ b/src/nomad_simulations/physical_property.py @@ -51,6 +51,9 @@ class PhysicalProperty(ArchiveSection): string identifiers and quantities for referencing and establishing the character of a physical property. """ + # TODO add `errors` + # TODO add `smearing` + name = Quantity( type=str, description=""" diff --git a/src/nomad_simulations/properties/__init__.py b/src/nomad_simulations/properties/__init__.py index 29011f91..c84a0fc8 100644 --- a/src/nomad_simulations/properties/__init__.py +++ b/src/nomad_simulations/properties/__init__.py @@ -16,4 +16,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from .energies import FermiLevel, ChemicalPotential from .band_gap import ElectronicBandGap +from .spectral_profile import ( + SpectralProfile, + DOSProfile, + ElectronicDensityOfStates, + XASSpectra, +) diff --git a/src/nomad_simulations/properties/band_gap.py b/src/nomad_simulations/properties/band_gap.py index e1a424b7..2949a7e6 100644 --- a/src/nomad_simulations/properties/band_gap.py +++ b/src/nomad_simulations/properties/band_gap.py @@ -31,6 +31,8 @@ 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( diff --git a/src/nomad_simulations/properties/energies.py b/src/nomad_simulations/properties/energies.py new file mode 100644 index 00000000..8d8209a3 --- /dev/null +++ b/src/nomad_simulations/properties/energies.py @@ -0,0 +1,79 @@ +# +# 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 ..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) diff --git a/src/nomad_simulations/properties/spectral_profile.py b/src/nomad_simulations/properties/spectral_profile.py new file mode 100644 index 00000000..943c3c6a --- /dev/null +++ b/src/nomad_simulations/properties/spectral_profile.py @@ -0,0 +1,604 @@ +# +# 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 + +from ..utils import get_sibling_section +from ..physical_property import PhysicalProperty +from ..variables import Energy2 as Energy +from ..atoms_state import AtomsState, OrbitalsState +from .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 _get_energy_points(self, logger: BoundLogger) -> Optional[pint.Quantity]: + """ + Gets the `grid_points` of the `Energy` variable if the required `Energy` variable is present in the `variables`. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[pint.Quantity]): The `grid_points` of the `Energy` variable. + """ + for var in self.variables: + if isinstance(var, Energy): + return var.grid_points + logger.error( + 'The required `Energy` variable is not present in the `variables`.' + ) + return None + + 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. + """ + # ! We need schema for `ElectronicEigenvalues` to store `highest_occupied` and `lowest_occupied` + # TODO improve and check this normalization after implementing `ElectronicEigenvalues` + + # 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_occupied` 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_occupied_energy = ( + eigenvalues.lowest_occupied if eigenvalues is not None else None + ) + # and set defaults for `highest_occupied_energy` and `lowest_occupied_energy` in `m_cache` + if highest_occupied_energy is not None: + self.m_cache['highest_occupied_energy'] = highest_occupied_energy + if lowest_occupied_energy is not None: + self.m_cache['lowest_occupied_energy'] = lowest_occupied_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_occupied_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_occupied_energy` stored + in `m_cache` from `resolve_energies_origin()`. If the difference of `highest_occupied_energy` and + `lowest_occupied_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_occupied_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 + + # 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=data[0].variables, + ) + 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 = self._get_energy_points(logger) + if energies is None: + return + + # 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 XASSpectra(SpectralProfile): + """ + X-ray Absorption Spectra (XAS). + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + xanes_spectra = SubSection( + sub_section=SpectralProfile.m_def, + description=""" + X-ray Absorption Near Edge Structure (XANES) spectra. + """, + repeats=False, + ) + + exafs_spectra = SubSection( + sub_section=SpectralProfile.m_def, + description=""" + Extended X-ray Absorption Fine Structure (EXAFS) spectra. + """, + 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 spectra 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_spectra is not None and self.exafs_spectra is not None: + # Concatenate XANE and EXAFS `Energy` grid points + xanes_energies = self.xanes_spectra._get_energy_points(logger) + exafs_energies = self.exafs_spectra._get_energy_points(logger) + if xanes_energies is None or exafs_energies is None: + logger.warning( + 'Could not extract the `Energy` grid points from XANES or EXAFS.' + ) + return + 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(grid_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_spectra.value, self.exafs_spectra.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 `XASSpectra.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) diff --git a/tests/conftest.py b/tests/conftest.py index f0dd8e07..e22b3c5b 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -18,12 +18,24 @@ import os import pytest +from typing import List, Optional from nomad.units import ureg -from nomad_simulations.outputs import Outputs, SCFOutputs +from . import logger + +from nomad_simulations import Simulation +from nomad_simulations.model_system import ModelSystem, AtomicCell +from nomad_simulations.atoms_state import AtomsState, OrbitalsState +from nomad_simulations.model_method import ModelMethod from nomad_simulations.numerical_settings import SelfConsistency -from nomad_simulations.properties import ElectronicBandGap +from nomad_simulations.outputs import Outputs, SCFOutputs +from nomad_simulations.variables import Energy2 as Energy +from nomad_simulations.properties import ( + ElectronicBandGap, + DOSProfile, + ElectronicDensityOfStates, +) if os.getenv('_PYTEST_RAISE', '0') != '0': @@ -36,7 +48,66 @@ def pytest_internalerror(excinfo): raise excinfo.value -def get_scf_electronic_band_gap_template(threshold_change: float = 1e-3) -> SCFOutputs: +def generate_simulation( + model_system: Optional[ModelSystem] = None, + model_method: Optional[ModelMethod] = None, + outputs: Optional[Outputs] = None, +) -> Simulation: + """ + Generate a `Simulation` section with the main sub-sections, `ModelSystem`, `ModelMethod`, and `Outputs`. If `ModelSystem` + and `Outputs` are set, then it adds `ModelSystem` as a reference in `Outputs`. + """ + simulation = Simulation() + if model_method is not None: + simulation.model_method.append(model_method) + if model_system is not None: + simulation.model_system.append(model_system) + if outputs is not None: + simulation.outputs.append(outputs) + outputs.model_system_ref = model_system + return simulation + + +def generate_model_system( + type: str = 'original', + positions: List[List[float]] = [[0, 0, 0], [0.5, 0.5, 0.5]], + chemical_symbols: List[str] = ['Ga', 'As'], + orbitals_symbols: List[List[str]] = [['s'], ['px', 'py']], +) -> Optional[ModelSystem]: + """ + Generate a `ModelSystem` section with the given parameters. + """ + if len(chemical_symbols) != len(orbitals_symbols): + return None + + model_system = ModelSystem() + atomic_cell = AtomicCell(type=type, positions=positions * ureg.meter) + model_system.cell.append(atomic_cell) + + # Add atoms_state to the model_system + atoms_state = [] + for element, orbitals in zip(chemical_symbols, orbitals_symbols): + orbitals_state = [] + for orbital in orbitals: + orbitals_state.append( + OrbitalsState( + l_quantum_symbol=orbital[0], ml_quantum_symbol=orbital[1:] + ) + ) # TODO add this split setter as part of the `OrbitalsState` methods + atom_state = AtomsState(chemical_symbol=element, orbitals_state=orbitals_state) + # and obtain the atomic number for each AtomsState + atom_state.resolve_chemical_symbol_and_number(logger) + atoms_state.append(atom_state) + atomic_cell.atoms_state = atoms_state + return model_system + + +def generate_scf_electronic_band_gap_template( + threshold_change: float = 1e-3, +) -> SCFOutputs: + """ + Generate a `SCFOutputs` section with a template for the electronic_band_gap property. + """ scf_outputs = SCFOutputs() # Define a list of scf_steps with values of the total energy like [1, 1.1, 1.11, 1.111, etc], # such that the difference between one step and the next one decreases a factor of 10. @@ -57,6 +128,56 @@ def get_scf_electronic_band_gap_template(threshold_change: float = 1e-3) -> SCFO return scf_outputs +def generate_simulation_electronic_dos( + energy_points: List[int] = [-3, -2, -1, 0, 1, 2, 3], +) -> Simulation: + """ + Generate a `Simulation` section with an `ElectronicDensityOfStates` section under `Outputs`. It uses + the template of the model_system created with the `generate_model_system` function. + """ + # Create the `Simulation` section to make refs work + model_system = generate_model_system() + outputs = Outputs() + simulation = generate_simulation(model_system=model_system, outputs=outputs) + + # Populating the `ElectronicDensityOfStates` section + variables_energy = [Energy(grid_points=energy_points * ureg.joule)] + electronic_dos = ElectronicDensityOfStates(variables=variables_energy) + outputs.electronic_dos.append(electronic_dos) + # electronic_dos.value = total_dos * ureg('1/joule') + orbital_s_Ga_pdos = DOSProfile( + variables=variables_energy, + entity_ref=model_system.cell[0].atoms_state[0].orbitals_state[0], + ) + orbital_px_As_pdos = DOSProfile( + variables=variables_energy, + entity_ref=model_system.cell[0].atoms_state[1].orbitals_state[0], + ) + orbital_py_As_pdos = DOSProfile( + variables=variables_energy, + entity_ref=model_system.cell[0].atoms_state[1].orbitals_state[1], + ) + orbital_s_Ga_pdos.value = [0.2, 0.5, 0, 0, 0, 0.0, 0.0] * ureg('1/joule') + orbital_px_As_pdos.value = [1.0, 0.2, 0, 0, 0, 0.3, 0.0] * ureg('1/joule') + orbital_py_As_pdos.value = [0.3, 0.5, 0, 0, 0, 0.5, 1.3] * ureg('1/joule') + electronic_dos.projected_dos = [ + orbital_s_Ga_pdos, + orbital_px_As_pdos, + orbital_py_As_pdos, + ] + return simulation + + +@pytest.fixture(scope='session') +def model_system() -> ModelSystem: + return generate_model_system() + + @pytest.fixture(scope='session') def scf_electronic_band_gap() -> SCFOutputs: - return get_scf_electronic_band_gap_template() + return generate_scf_electronic_band_gap_template() + + +@pytest.fixture(scope='session') +def simulation_electronic_dos() -> Simulation: + return generate_simulation_electronic_dos() diff --git a/tests/test_energies.py b/tests/test_energies.py new file mode 100644 index 00000000..0c6f9722 --- /dev/null +++ b/tests/test_energies.py @@ -0,0 +1,60 @@ +# +# 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 pytest +import numpy as np +from typing import Optional, List, Union + +from . import logger + +from nomad_simulations.properties import FermiLevel, ChemicalPotential + + +class TestFermiLevel: + """ + Test the `FermiLevel` class defined in `properties/energies.py`. + """ + + # ! Include this initial `test_default_quantities` method when testing your PhysicalProperty classes + def test_default_quantities(self): + """ + Test the default quantities assigned when creating an instance of the `FermiLevel` class. + """ + fermi_level = FermiLevel() + assert fermi_level.iri == 'http://fairmat-nfdi.eu/taxonomy/FermiLevel' + assert fermi_level.name == 'FermiLevel' + assert fermi_level.rank == [] + + +class TestChemicalPotential: + """ + Test the `ChemicalPotential` class defined in `properties/energies.py`. + """ + + # ! Include this initial `test_default_quantities` method when testing your PhysicalProperty classes + def test_default_quantities(self): + """ + Test the default quantities assigned when creating an instance of the `ChemicalPotential` class. + """ + chemical_potential = ChemicalPotential() + assert ( + chemical_potential.iri + == 'http://fairmat-nfdi.eu/taxonomy/ChemicalPotential' + ) + assert chemical_potential.name == 'ChemicalPotential' + assert chemical_potential.rank == [] diff --git a/tests/test_outputs.py b/tests/test_outputs.py index 5e5befa7..756f6b8b 100644 --- a/tests/test_outputs.py +++ b/tests/test_outputs.py @@ -17,31 +17,14 @@ # import pytest -import numpy as np from . import logger -from .conftest import get_scf_electronic_band_gap_template +from .conftest import generate_scf_electronic_band_gap_template from nomad.units import ureg -from nomad.metainfo import Quantity -from nomad_simulations.physical_property import PhysicalProperty from nomad_simulations.outputs import Outputs, ElectronicBandGap -class TotalEnergy(PhysicalProperty): - """Physical property class defined for testing purposes.""" - - rank = [] - variables = [] - value = Quantity( - type=np.float64, - unit='joule', - description=""" - The total energy of the system. - """, - ) - - class TestOutputs: """ Test the `Outputs` class defined in `outputs.py`. @@ -55,7 +38,7 @@ def test_is_scf_converged(self, threshold_change: float, result: bool): """ Test the `resolve_is_scf_converged` method. """ - scf_outputs = get_scf_electronic_band_gap_template( + scf_outputs = generate_scf_electronic_band_gap_template( threshold_change=threshold_change ) is_scf_converged = scf_outputs.resolve_is_scf_converged( @@ -99,14 +82,9 @@ def test_normalize(self, threshold_change: float, result: bool): """ Test the `normalize` method. """ - scf_outputs = get_scf_electronic_band_gap_template( + scf_outputs = generate_scf_electronic_band_gap_template( threshold_change=threshold_change ) - # Add a non-SCF calculated PhysicalProperty - scf_outputs.custom_physical_property = [ - TotalEnergy(name='TotalEnergy', value=1 * ureg.joule) - ] scf_outputs.normalize(None, logger) assert scf_outputs.electronic_band_gap[0].is_scf_converged == result - assert scf_outputs.custom_physical_property[0].is_scf_converged is None diff --git a/tests/test_spectral_profile.py b/tests/test_spectral_profile.py new file mode 100644 index 00000000..c28b0989 --- /dev/null +++ b/tests/test_spectral_profile.py @@ -0,0 +1,303 @@ +# +# 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 pytest +import numpy as np +from typing import List, Optional + +from . import logger + +from nomad.units import ureg + +from nomad_simulations import Simulation +from nomad_simulations.model_system import ModelSystem, AtomicCell +from nomad_simulations.atoms_state import AtomsState +from nomad_simulations.outputs import Outputs +from nomad_simulations.properties import ( + SpectralProfile, + ElectronicDensityOfStates, + XASSpectra, +) +from nomad_simulations.variables import Temperature, Energy2 as Energy + + +class TestSpectralProfile: + """ + Test the `SpectralProfile` class defined in `properties/spectral_profile.py`. + """ + + def test_is_valid_spectral_profile(self): + """ + Test the `is_valid_spectral_profile` method. + """ + spectral_profile = SpectralProfile( + variables=[Energy(grid_points=[-1, 0, 1] * ureg.joule)] + ) + spectral_profile.value = [1.5, 0, 0.8] + assert spectral_profile.is_valid_spectral_profile() + spectral_profile.value = [2, 0, -4] + assert not spectral_profile.is_valid_spectral_profile() + + +class TestElectronicDensityOfStates: + """ + Test the `ElectronicDensityOfStates` class defined in `properties/spectral_profile.py`. + """ + + # ! Include this initial `test_default_quantities` method when testing your PhysicalProperty classes + def test_default_quantities(self): + """ + Test the default quantities assigned when creating an instance of the `ElectronicDensityOfStates` class. + """ + electronic_dos = ElectronicDensityOfStates() + assert ( + electronic_dos.iri + == 'http://fairmat-nfdi.eu/taxonomy/ElectronicDensityOfStates' + ) + assert electronic_dos.name == 'ElectronicDensityOfStates' + assert electronic_dos.rank == [] + + def test_get_energy_points(self): + """ + Test the `_get_energy_points` private method. We test here that the function indeed returns the points from a `Energy` variable. + """ + electronic_dos = ElectronicDensityOfStates() + electronic_dos.variables = [ + Temperature(grid_points=list(range(-3, 4)) * ureg.kelvin) + ] + assert electronic_dos._get_energy_points(logger) is None + electronic_dos.variables.append( + Energy(grid_points=list(range(-3, 4)) * ureg.joule) + ) + energies = electronic_dos._get_energy_points(logger) + assert (energies.magnitude == np.array([-3, -2, -1, 0, 1, 2, 3])).all() + + def test_resolve_energies_origin(self): + """ + Test the `resolve_energies_origin` method. + """ + # ! add test when `ElectronicEigenvalues` is implemented + pass + + def test_resolve_normalization_factor(self, simulation_electronic_dos: Simulation): + """ + Test the `resolve_normalization_factor` method. + """ + simulation = Simulation() + outputs = Outputs() + # We only used the `simulation_electronic_dos` fixture to get the `ElectronicDensityOfStates` to test missing refs + electronic_dos = simulation_electronic_dos.outputs[0].electronic_dos[0] + electronic_dos.energies_origin = 0.5 * ureg.joule + outputs.electronic_dos.append(electronic_dos) + simulation.outputs.append(outputs) + + # No `model_system_ref` + assert electronic_dos.resolve_normalization_factor(logger) is None + + # No `model_system_ref.cell` + model_system = ModelSystem() + simulation.model_system.append(model_system) + outputs.model_system_ref = simulation.model_system[0] + assert electronic_dos.resolve_normalization_factor(logger) is None + + # No `model_system_ref.cell.atoms_state` + atomic_cell = AtomicCell( + type='original', positions=[[0, 0, 0], [0.5, 0.5, 0.5]] * ureg.meter + ) + model_system.cell.append(atomic_cell) + assert electronic_dos.resolve_normalization_factor(logger) is None + + # Adding the required `model_system_ref` sections and quantities + atoms_state = [ + AtomsState(chemical_symbol='Ga'), + AtomsState(chemical_symbol='As'), + ] + for atom in atoms_state: + atom.resolve_chemical_symbol_and_number(logger) + atomic_cell.atoms_state = atoms_state + # Non spin-polarized + normalization_factor = electronic_dos.resolve_normalization_factor(logger) + assert np.isclose(normalization_factor, 0.015625) + # Spin-polarized + electronic_dos.spin_channel = 0 + normalization_factor_spin_polarized = ( + electronic_dos.resolve_normalization_factor(logger) + ) + assert np.isclose( + normalization_factor_spin_polarized, 0.5 * normalization_factor + ) + + def test_extract_band_gap(self): + """ + Test the `extract_band_gap` method. + """ + # ! add test when `ElectronicEigenvalues` is implemented + pass + + def test_resolve_pdos_name(self, simulation_electronic_dos: Simulation): + """ + Test the `resolve_pdos_name` method. + """ + # Get projected DOSProfile from the simulation fixture + projected_dos = ( + simulation_electronic_dos.outputs[0].electronic_dos[0].projected_dos + ) + assert len(projected_dos) == 3 + pdos_names = ['orbital s Ga', 'orbital px As', 'orbital py As'] + for i, pdos in enumerate(projected_dos): + name = pdos.resolve_pdos_name(logger) + assert name == pdos_names[i] + + def test_extract_projected_dos(self, simulation_electronic_dos: Simulation): + """ + Test the `extract_projected_dos` method. + """ + # Get Outputs and ElectronicDensityOfStates from the simulation fixture + outputs = simulation_electronic_dos.outputs[0] + electronic_dos = outputs.electronic_dos[0] + + # Initial tests for the passed `projected_dos` (only orbital PDOS) + assert len(electronic_dos.projected_dos) == 3 # only orbital projected DOS + orbital_projected = electronic_dos.extract_projected_dos('orbital', logger) + atom_projected = electronic_dos.extract_projected_dos('atom', logger) + assert len(orbital_projected) == 3 and len(atom_projected) == 0 + orbital_projected_names = [orb_pdos.name for orb_pdos in orbital_projected] + assert orbital_projected_names == [ + 'orbital s Ga', + 'orbital px As', + 'orbital py As', + ] + assert ( + orbital_projected[0].entity_ref + == outputs.model_system_ref.cell[0].atoms_state[0].orbitals_state[0] + ) # orbital `s` in `Ga` atom + assert ( + orbital_projected[1].entity_ref + == outputs.model_system_ref.cell[0].atoms_state[1].orbitals_state[0] + ) # orbital `px` in `As` atom + assert ( + orbital_projected[2].entity_ref + == outputs.model_system_ref.cell[0].atoms_state[1].orbitals_state[1] + ) # orbital `py` in `As` atom + + orbital_projected = electronic_dos.extract_projected_dos('orbital', logger) + atom_projected = electronic_dos.extract_projected_dos('atom', logger) + assert len(orbital_projected) == 3 and len(atom_projected) == 0 + + @pytest.mark.parametrize( + 'value, result', + [ + (None, [1.5, 1.2, 0, 0, 0, 0.8, 1.3]), + ([30.5, 1.2, 0, 0, 0, 0.8, 1.3], [30.5, 1.2, 0, 0, 0, 0.8, 1.3]), + ], + ) + def test_generate_from_pdos( + self, + simulation_electronic_dos: Simulation, + value: Optional[List[float]], + result: List[float], + ): + """ + Test the `generate_from_projected_dos` method. + """ + # Get Outputs and ElectronicDensityOfStates from the simulation fixture + outputs = simulation_electronic_dos.outputs[0] + electronic_dos = outputs.electronic_dos[0] + + # Add `value` + if value is not None: + electronic_dos.value = value * ureg('1/joule') + + val = electronic_dos.generate_from_projected_dos(logger) + assert (val.magnitude == result).all() + + # Testing orbital + atom projected DOS (3 orbitals + 2 atoms PDOS) + assert len(electronic_dos.projected_dos) == 5 + orbital_projected = electronic_dos.extract_projected_dos('orbital', logger) + atom_projected = electronic_dos.extract_projected_dos('atom', logger) + assert len(orbital_projected) == 3 and len(atom_projected) == 2 + atom_projected_names = [atom_pdos.name for atom_pdos in atom_projected] + assert atom_projected_names == ['atom Ga', 'atom As'] + assert ( + atom_projected[0].entity_ref + == outputs.model_system_ref.cell[0].atoms_state[0] + ) # `Ga` atom + assert ( + atom_projected[1].entity_ref + == outputs.model_system_ref.cell[0].atoms_state[1] + ) # `As` atom + + def test_normalize(self): + """ + Test the `normalize` method. + """ + # ! add test when `ElectronicEigenvalues` is implemented + pass + + +class TestXASSpectra: + """ + Test the `XASSpectra` class defined in `properties/spectral_profile.py`. + """ + + # ! Include this initial `test_default_quantities` method when testing your PhysicalProperty classes + def test_default_quantities(self): + """ + Test the default quantities assigned when creating an instance of the `XASSpectra` class. + """ + xas_spectra = XASSpectra() + assert xas_spectra.iri is None # Add iri when available + assert xas_spectra.name == 'XASSpectra' + assert xas_spectra.rank == [] + + @pytest.mark.parametrize( + 'xanes_energies, exafs_energies, xas_values', + [ + (None, None, None), + ([0, 1, 2], None, None), + (None, [3, 4, 5], None), + ([0, 1, 2], [3, 4, 5], [0.5, 0.1, 0.3, 0.2, 0.4, 0.6]), + ([0, 1, 4], [3, 4, 5], None), + ([0, 1, 2], [0, 4, 5], None), + ], + ) + def test_generate_from_contributions( + self, + xanes_energies: Optional[List[float]], + exafs_energies: Optional[List[float]], + xas_values: Optional[List[float]], + ): + """ + Test the `generate_from_contributions` method. + """ + xas_spectra = XASSpectra() + if xanes_energies is not None: + xanes_spectra = SpectralProfile() + xanes_spectra.variables = [Energy(grid_points=xanes_energies * ureg.joule)] + xanes_spectra.value = [0.5, 0.1, 0.3] + xas_spectra.xanes_spectra = xanes_spectra + if exafs_energies is not None: + exafs_spectra = SpectralProfile() + exafs_spectra.variables = [Energy(grid_points=exafs_energies * ureg.joule)] + exafs_spectra.value = [0.2, 0.4, 0.6] + xas_spectra.exafs_spectra = exafs_spectra + xas_spectra.generate_from_contributions(logger) + if xas_spectra.value is not None: + assert (xas_spectra.value == xas_values).all() + else: + assert xas_spectra.value == xas_values