From c1459240fee0594605ea1bbb9d2e0472c14ed76b Mon Sep 17 00:00:00 2001 From: JosePizarro3 Date: Mon, 27 May 2024 16:57:23 +0200 Subject: [PATCH] Squashing stuff to avoid issues when rebasing... --- src/nomad_simulations/numerical_settings.py | 94 +++++++++++++ src/nomad_simulations/outputs.py | 20 ++- src/nomad_simulations/properties/__init__.py | 4 +- .../properties/permittivity.py | 118 ++++++++++++++++ .../properties/spectral_profile.py | 93 +++++++----- src/nomad_simulations/utils/__init__.py | 7 +- src/nomad_simulations/utils/utils.py | 25 +++- src/nomad_simulations/variables.py | 22 +++ tests/conftest.py | 6 +- tests/test_numerical_settings.py | 2 - tests/test_outputs.py | 16 +-- tests/test_permittivity.py | 133 ++++++++++++++++++ tests/test_spectral_profile.py | 74 +++++----- tests/test_utils.py | 34 ++++- 14 files changed, 551 insertions(+), 97 deletions(-) create mode 100644 src/nomad_simulations/properties/permittivity.py create mode 100644 tests/test_permittivity.py diff --git a/src/nomad_simulations/numerical_settings.py b/src/nomad_simulations/numerical_settings.py index 6cfabcb5..5954d46a 100644 --- a/src/nomad_simulations/numerical_settings.py +++ b/src/nomad_simulations/numerical_settings.py @@ -462,6 +462,99 @@ def resolve_k_line_density( return k_line_density return None + def resolve_high_symmetry_points( + self, + model_systems: List[ModelSystem], + logger: BoundLogger, + eps: float = 3e-3, + ) -> Optional[dict]: + """ + Resolves the `high_symmetry_points` of the `KMesh` from the list of `ModelSystem`. This method + relies on using the `ModelSystem` information in the sub-sections `Symmetry` and `AtomicCell`, and uses + the ASE package to extract the special (high symmetry) points information. + + Args: + model_systems (List[ModelSystem]): The list of `ModelSystem` sections. + logger (BoundLogger): The logger to log messages. + eps (float, optional): Tolerance factor to define the `lattice` ASE object. Defaults to 3e-3. + + Returns: + (Optional[dict]): The resolved `high_symmetry_points` of the `KMesh`. + """ + # Extracting `bravais_lattice` from `ModelSystem.symmetry` section and `ASE.cell` from `ModelSystem.cell` + lattice = None + for model_system in model_systems: + # General checks to proceed with normalization + if is_not_representative(model_system, logger): + continue + if model_system.symmetry is None: + logger.warning('Could not find `ModelSystem.symmetry`.') + continue + bravais_lattice = [symm.bravais_lattice for symm in model_system.symmetry] + if len(bravais_lattice) != 1: + logger.warning( + 'Could not uniquely determine `bravais_lattice` from `ModelSystem.symmetry`.' + ) + continue + bravais_lattice = bravais_lattice[0] + + if model_system.cell is None: + logger.warning('Could not find `ModelSystem.cell`.') + continue + prim_atomic_cell = None + for atomic_cell in model_system.cell: + if atomic_cell.type == 'primitive': + prim_atomic_cell = atomic_cell + break + if prim_atomic_cell is None: + logger.warning( + 'Could not find the primitive `AtomicCell` under `ModelSystem.cell`.' + ) + continue + # function defined in AtomicCell + atoms = prim_atomic_cell.to_ase_atoms(logger) + cell = atoms.get_cell() + lattice = cell.get_bravais_lattice(eps) + break # only cover the first representative `ModelSystem` + + # Checking if `bravais_lattice` and `lattice` are defined + if lattice is None: + logger.warning( + 'Could not resolve `bravais_lattice` and `lattice` ASE object from the `ModelSystem`.' + ) + return None + + # Non-conventional ordering testing for certain lattices: + if bravais_lattice in ['oP', 'oF', 'oI', 'oS']: + a, b, c = lattice.a, lattice.b, lattice.c + assert a < b + if bravais_lattice != 'oS': + assert b < c + elif bravais_lattice in ['mP', 'mS']: + a, b, c = lattice.a, lattice.b, lattice.c + alpha = lattice.alpha * np.pi / 180 + assert a <= c and b <= c # ordering of the conventional lattice + assert alpha < np.pi / 2 + + # Extracting the `high_symmetry_points` from the `lattice` object + special_points = lattice.get_special_points() + if special_points is None: + logger.warning( + 'Could not find `lattice.get_special_points()` from the ASE package.' + ) + return None + high_symmetry_points = {} + for key, value in lattice.get_special_points().items(): + if key == 'G': + key = 'Gamma' + if bravais_lattice == 'tI': + if key == 'S': + key = 'Sigma' + elif key == 'S1': + key = 'Sigma1' + high_symmetry_points[key] = list(value) + return high_symmetry_points + def normalize(self, archive, logger) -> None: super().normalize(archive, logger) @@ -679,6 +772,7 @@ def resolve_points( 'The `reciprocal_lattice_vectors` are not passed as an input.' ) return None + # Check if `points_norm` is a list and convert it to a numpy array if isinstance(points_norm, list): points_norm = np.array(points_norm) diff --git a/src/nomad_simulations/outputs.py b/src/nomad_simulations/outputs.py index ac1ca674..46436a6f 100644 --- a/src/nomad_simulations/outputs.py +++ b/src/nomad_simulations/outputs.py @@ -33,7 +33,9 @@ HoppingMatrix, ElectronicBandGap, ElectronicDensityOfStates, - XASSpectra, + AbsorptionSpectrum, + XASSpectrum, + Permittivity, ) @@ -63,23 +65,27 @@ class Outputs(ArchiveSection): # List of properties # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # - fermi_level = SubSection(sub_section=FermiLevel.m_def, repeats=True) + fermi_levels = SubSection(sub_section=FermiLevel.m_def, repeats=True) - chemical_potential = SubSection(sub_section=ChemicalPotential.m_def, repeats=True) + chemical_potentials = SubSection(sub_section=ChemicalPotential.m_def, repeats=True) - crystal_field_splitting = SubSection( + crystal_field_splittings = SubSection( sub_section=CrystalFieldSplitting.m_def, repeats=True ) - hopping_matrix = SubSection(sub_section=HoppingMatrix.m_def, repeats=True) + hopping_matrices = SubSection(sub_section=HoppingMatrix.m_def, repeats=True) - electronic_band_gap = SubSection(sub_section=ElectronicBandGap.m_def, repeats=True) + electronic_band_gaps = 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) + permittivities = SubSection(sub_section=Permittivity.m_def, repeats=True) + + absorption_spectra = SubSection(sub_section=AbsorptionSpectrum.m_def, repeats=True) + + xas_spectra = SubSection(sub_section=XASSpectrum.m_def, repeats=True) # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # diff --git a/src/nomad_simulations/properties/__init__.py b/src/nomad_simulations/properties/__init__.py index e529ed38..889196cd 100644 --- a/src/nomad_simulations/properties/__init__.py +++ b/src/nomad_simulations/properties/__init__.py @@ -22,6 +22,8 @@ SpectralProfile, DOSProfile, ElectronicDensityOfStates, - XASSpectra, + AbsorptionSpectrum, + XASSpectrum, ) from .hopping_matrix import HoppingMatrix, CrystalFieldSplitting +from .permittivity import Permittivity diff --git a/src/nomad_simulations/properties/permittivity.py b/src/nomad_simulations/properties/permittivity.py new file mode 100644 index 00000000..f4e4511c --- /dev/null +++ b/src/nomad_simulations/properties/permittivity.py @@ -0,0 +1,118 @@ +# +# 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/src/nomad_simulations/properties/spectral_profile.py b/src/nomad_simulations/properties/spectral_profile.py index 4b57ba60..03bd60ac 100644 --- a/src/nomad_simulations/properties/spectral_profile.py +++ b/src/nomad_simulations/properties/spectral_profile.py @@ -22,9 +22,9 @@ import pint from nomad import config -from nomad.metainfo import Quantity, SubSection, Section, Context +from nomad.metainfo import Quantity, SubSection, Section, Context, MEnum -from ..utils import get_sibling_section +from ..utils import get_sibling_section, get_variables from ..physical_property import PhysicalProperty from ..variables import Energy2 as Energy from ..atoms_state import AtomsState, OrbitalsState @@ -49,24 +49,6 @@ def __init__( super().__init__(m_def, m_context, **kwargs) self.rank = [] - def _get_energy_points(self, logger: BoundLogger) -> Optional[pint.Quantity]: - """ - Gets the `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 `points` of the `Energy` variable. - """ - for var in self.variables: - if isinstance(var, Energy): - return var.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. @@ -450,6 +432,14 @@ def generate_from_projected_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) @@ -468,7 +458,7 @@ def generate_from_projected_dos( atom_dos = DOSProfile( name=f'atom {ref.chemical_symbol}', entity_ref=ref, - variables=data[0].variables, + variables=energies, ) orbital_values = [ dos.value.magnitude for dos in data @@ -493,9 +483,10 @@ 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: + 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( @@ -528,25 +519,49 @@ def normalize(self, archive, logger) -> None: self.value = value_from_pdos -class XASSpectra(SpectralProfile): +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 Spectra (XAS). + X-ray Absorption Spectrum (XAS). """ # ! implement `iri` and `rank` as part of `m_def = Section()` - xanes_spectra = SubSection( - sub_section=SpectralProfile.m_def, + xanes_spectrum = SubSection( + sub_section=AbsorptionSpectrum.m_def, description=""" - X-ray Absorption Near Edge Structure (XANES) spectra. + X-ray Absorption Near Edge Structure (XANES) spectrum. """, repeats=False, ) - exafs_spectra = SubSection( - sub_section=SpectralProfile.m_def, + exafs_spectrum = SubSection( + sub_section=AbsorptionSpectrum.m_def, description=""" - Extended X-ray Absorption Fine Structure (EXAFS) spectra. + Extended X-ray Absorption Fine Structure (EXAFS) spectrum. """, repeats=False, ) @@ -560,22 +575,24 @@ def __init__( 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 + 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_spectra is not None and self.exafs_spectra is not None: + if self.xanes_spectrum is not None and self.exafs_spectrum 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: + 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.' @@ -587,7 +604,7 @@ def generate_from_contributions(self, logger: BoundLogger) -> None: # 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] + [self.xanes_spectrum.value, self.exafs_spectrum.value] ) except ValueError: logger.warning( @@ -599,6 +616,6 @@ def normalize(self, archive, logger) -> None: 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.' + '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) diff --git a/src/nomad_simulations/utils/__init__.py b/src/nomad_simulations/utils/__init__.py index a0a83520..0c96821f 100644 --- a/src/nomad_simulations/utils/__init__.py +++ b/src/nomad_simulations/utils/__init__.py @@ -16,4 +16,9 @@ # See the License for the specific language governing permissions and # limitations under the License. -from .utils import get_sibling_section, RussellSaundersState, is_not_representative +from .utils import ( + get_sibling_section, + RussellSaundersState, + is_not_representative, + get_variables, +) diff --git a/src/nomad_simulations/utils/utils.py b/src/nomad_simulations/utils/utils.py index e2fe582f..e48a5118 100644 --- a/src/nomad_simulations/utils/utils.py +++ b/src/nomad_simulations/utils/utils.py @@ -18,7 +18,7 @@ # from math import factorial -from typing import Optional +from typing import Optional, List from structlog.stdlib import BoundLogger from nomad.datamodel.data import ArchiveSection @@ -129,3 +129,26 @@ def is_not_representative(model_system, logger: BoundLogger): logger.warning('The `ModelSystem` was not found to be representative.') return True return False + + +# cannot define typing with `Variables` due to circular import issue +def get_variables( + variables: Optional[List[ArchiveSection]], variable_cls: ArchiveSection +) -> List[ArchiveSection]: + """ + Get the list of variables which are of type `variable_cls` and appear under `variables`. + + Args: + variables (List[Variables]): The list of variables to check. + variable_cls (Variables): The class of the variables to get. + + Returns: + (List[Variables]): The list of variables which are of type `variable_cls`. + """ + if variables is None: + return [] + result = [] + for var in variables: + if isinstance(var, variable_cls): + result.append(var) + return result diff --git a/src/nomad_simulations/variables.py b/src/nomad_simulations/variables.py index 3a786776..166ea90e 100644 --- a/src/nomad_simulations/variables.py +++ b/src/nomad_simulations/variables.py @@ -157,6 +157,28 @@ def normalize(self, archive, logger) -> None: super().normalize(archive, logger) +class Frequency(Variables): + """ """ + + points = Quantity( + type=np.float64, + unit='joule', + shape=['n_points'], + description=""" + Points in which the frequency is discretized in joules. + """, + ) + + 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) + + class KMesh(Variables): """ K-point mesh over which the physical property is calculated. This is used to define `ElectronicEigenvalues(PhysicalProperty)` and diff --git a/tests/conftest.py b/tests/conftest.py index 39107094..2f62bb31 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -170,16 +170,16 @@ def generate_scf_electronic_band_gap_template( for i in range(1, n_scf_steps): value = 1 + sum([1 / (10**j) for j in range(1, i + 1)]) scf_step = Outputs( - electronic_band_gap=[ElectronicBandGap(value=value * ureg.joule)] + electronic_band_gaps=[ElectronicBandGap(value=value * ureg.joule)] ) scf_outputs.scf_steps.append(scf_step) # Add a SCF calculated PhysicalProperty - scf_outputs.electronic_band_gap.append(ElectronicBandGap(value=value * ureg.joule)) + scf_outputs.electronic_band_gaps.append(ElectronicBandGap(value=value * ureg.joule)) # and a `SelfConsistency` ref section scf_params = SelfConsistency( threshold_change=threshold_change, threshold_change_unit='joule' ) - scf_outputs.electronic_band_gap[0].self_consistency_ref = scf_params + scf_outputs.electronic_band_gaps[0].self_consistency_ref = scf_params return scf_outputs diff --git a/tests/test_numerical_settings.py b/tests/test_numerical_settings.py index 3b780da7..11d94ba4 100644 --- a/tests/test_numerical_settings.py +++ b/tests/test_numerical_settings.py @@ -252,7 +252,6 @@ def test_resolve_k_line_density( reciprocal_lattice_vectors = k_space.reciprocal_lattice_vectors k_mesh = k_space.k_mesh[0] model_systems = simulation.model_system - # Applying method `get_k_line_density` get_k_line_density_value = k_mesh.get_k_line_density( reciprocal_lattice_vectors=reciprocal_lattice_vectors, logger=logger @@ -264,7 +263,6 @@ def test_resolve_k_line_density( ) else: assert get_k_line_density_value == result_get_k_line_density - # Applying method `resolve_k_line_density` k_line_density = k_mesh.resolve_k_line_density( model_systems=model_systems, diff --git a/tests/test_outputs.py b/tests/test_outputs.py index aac05ced..d4efaaa2 100644 --- a/tests/test_outputs.py +++ b/tests/test_outputs.py @@ -44,9 +44,9 @@ def test_is_scf_converged(self, threshold_change: float, result: bool): threshold_change=threshold_change ) is_scf_converged = scf_outputs.resolve_is_scf_converged( - property_name='electronic_band_gap', + property_name='electronic_band_gaps', i_property=0, - phys_property=scf_outputs.electronic_band_gap[0], + phys_property=scf_outputs.electronic_band_gaps[0], logger=logger, ) assert is_scf_converged == result @@ -60,18 +60,18 @@ def test_extract_spin_polarized_properties(self): # No spin-polarized band gap band_gap_non_spin_polarized = ElectronicBandGap(variables=[]) band_gap_non_spin_polarized.value = 2.0 * ureg.joule - outputs.electronic_band_gap.append(band_gap_non_spin_polarized) - band_gaps = outputs.extract_spin_polarized_property('electronic_band_gap') + outputs.electronic_band_gaps.append(band_gap_non_spin_polarized) + band_gaps = outputs.extract_spin_polarized_property('electronic_band_gaps') assert band_gaps == [] # Spin-polarized band gaps band_gap_spin_1 = ElectronicBandGap(variables=[], spin_channel=0) band_gap_spin_1.value = 1.0 * ureg.joule - outputs.electronic_band_gap.append(band_gap_spin_1) + outputs.electronic_band_gaps.append(band_gap_spin_1) band_gap_spin_2 = ElectronicBandGap(variables=[], spin_channel=1) band_gap_spin_2.value = 1.5 * ureg.joule - outputs.electronic_band_gap.append(band_gap_spin_2) - band_gaps = outputs.extract_spin_polarized_property('electronic_band_gap') + outputs.electronic_band_gaps.append(band_gap_spin_2) + band_gaps = outputs.extract_spin_polarized_property('electronic_band_gaps') assert len(band_gaps) == 2 assert band_gaps[0].value.magnitude == 1.0 assert band_gaps[1].value.magnitude == 1.5 @@ -89,4 +89,4 @@ def test_normalize(self, threshold_change: float, result: bool): ) scf_outputs.normalize(EntryArchive(), logger) - assert scf_outputs.electronic_band_gap[0].is_scf_converged == result + assert scf_outputs.electronic_band_gaps[0].is_scf_converged == result diff --git a/tests/test_permittivity.py b/tests/test_permittivity.py new file mode 100644 index 00000000..23743588 --- /dev/null +++ b/tests/test_permittivity.py @@ -0,0 +1,133 @@ +# +# 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 nomad.units import ureg +from nomad.datamodel import EntryArchive + +from . import logger + +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 Permittivity +from nomad_simulations.variables import Variables, KMesh, Frequency + + +class TestPermittivity: + """ + Test the `Permittivity` class defined in `properties/permittivity.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 `Permittivity` class. + """ + permittivity = Permittivity() + assert permittivity.iri == 'http://fairmat-nfdi.eu/taxonomy/Permittivity' + assert permittivity.name == 'Permittivity' + assert permittivity.rank == [3, 3] + + @pytest.mark.parametrize( + 'variables, result', + [ + (None, 'static'), + ([], 'static'), + ([KMesh()], 'static'), + ([KMesh(), Frequency()], 'dynamic'), + ], + ) + def test_resolve_type(self, variables: Optional[List[Variables]], result: str): + """ + Test the `resolve_type` method. + """ + permittivity = Permittivity() + if variables is not None: + permittivity.variables = [var for var in variables] + assert permittivity.resolve_type() == result + + @pytest.mark.parametrize( + 'variables, value, result', + [ + # Empty case + (None, None, None), + # No `variables` + ([], np.eye(3) * (1 + 1j), None), + # If `variables` contain `KMesh`, we cannot extract absorption spectra + ( + [KMesh(points=[[1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4]])], + np.array([np.eye(3) * k_point * (1 + 1j) for k_point in range(1, 5)]), + None, + ), + # Even if `variables` contain `Frequency`, we cannot extract absorption spectra if `value` depends on `KMesh` + ( + [ + Frequency(points=[0, 1, 2, 3, 4]), + KMesh(points=[[1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4]]), + ], + np.array( + [ + [ + np.eye(3) * k_point * (1 + 1j) + + np.eye(3) * freq_point * 0.5j + for k_point in range(1, 5) + ] + for freq_point in range(5) + ] + ), + None, + ), + # Valid case: `value` does not depend on `KMesh` and we can extract absorption spectra + ( + [ + Frequency(points=[0, 1, 2, 3, 4]), + ], + np.array([np.eye(3) * freq_point * 0.5j for freq_point in range(5)]), + [0.0, 0.5, 1.0, 1.5, 2.0], + ), + ], + ) + def test_extract_absorption_spectra( + self, + variables: Optional[List[Variables]], + value: Optional[np.ndarray], + result: Optional[List[float]], + ): + """ + Test the `extract_absorption_spectra` method. The `result` in the last valid case corresponds to the imaginary part of + the diagonal of the `Permittivity.value` for each frequency point. + """ + permittivity = Permittivity() + if variables is not None: + permittivity.variables = [var for var in variables] + permittivity.value = value + absorption_spectra = permittivity.extract_absorption_spectra(logger) + if absorption_spectra is not None: + assert len(absorption_spectra) == 3 + spectrum = absorption_spectra[1] + assert spectrum.rank == [] + assert spectrum.axis == 'yy' + assert len(spectrum.value) == len(variables[0].points) + assert np.allclose(spectrum.value, result) + else: + assert absorption_spectra == result diff --git a/tests/test_spectral_profile.py b/tests/test_spectral_profile.py index 4316ecf7..14a55034 100644 --- a/tests/test_spectral_profile.py +++ b/tests/test_spectral_profile.py @@ -32,9 +32,10 @@ from nomad_simulations.properties import ( SpectralProfile, ElectronicDensityOfStates, - XASSpectra, + AbsorptionSpectrum, + XASSpectrum, ) -from nomad_simulations.variables import Temperature, Energy2 as Energy +from nomad_simulations.variables import Energy2 as Energy class TestSpectralProfile: @@ -73,19 +74,6 @@ def test_default_quantities(self): 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(points=list(range(-3, 4)) * ureg.kelvin) - ] - assert electronic_dos._get_energy_points(logger) is None - electronic_dos.variables.append(Energy(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. @@ -249,20 +237,36 @@ def test_normalize(self): pass -class TestXASSpectra: +class TestAbsorptionSpectrum: + """ + Test the `AbsorptionSpectrum` 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 `AbsorptionSpectrum` class. + """ + absorption_spectrum = AbsorptionSpectrum() + assert absorption_spectrum.iri is None # Add iri when available + assert absorption_spectrum.name == 'AbsorptionSpectrum' + assert absorption_spectrum.rank == [] + + +class TestXASSpectrum: """ - Test the `XASSpectra` class defined in `properties/spectral_profile.py`. + Test the `XASSpectrum` 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. + Test the default quantities assigned when creating an instance of the `XASSpectrum` class. """ - xas_spectra = XASSpectra() - assert xas_spectra.iri is None # Add iri when available - assert xas_spectra.name == 'XASSpectra' - assert xas_spectra.rank == [] + xas_spectrum = XASSpectrum() + assert xas_spectrum.iri is None # Add iri when available + assert xas_spectrum.name == 'XASSpectrum' + assert xas_spectrum.rank == [] @pytest.mark.parametrize( 'xanes_energies, exafs_energies, xas_values', @@ -284,19 +288,19 @@ def test_generate_from_contributions( """ Test the `generate_from_contributions` method. """ - xas_spectra = XASSpectra() + xas_spectrum = XASSpectrum() if xanes_energies is not None: - xanes_spectra = SpectralProfile() - xanes_spectra.variables = [Energy(points=xanes_energies * ureg.joule)] - xanes_spectra.value = [0.5, 0.1, 0.3] - xas_spectra.xanes_spectra = xanes_spectra + xanes_spectrum = AbsorptionSpectrum() + xanes_spectrum.variables = [Energy(points=xanes_energies * ureg.joule)] + xanes_spectrum.value = [0.5, 0.1, 0.3] + xas_spectrum.xanes_spectrum = xanes_spectrum if exafs_energies is not None: - exafs_spectra = SpectralProfile() - exafs_spectra.variables = [Energy(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() + exafs_spectrum = AbsorptionSpectrum() + exafs_spectrum.variables = [Energy(points=exafs_energies * ureg.joule)] + exafs_spectrum.value = [0.2, 0.4, 0.6] + xas_spectrum.exafs_spectrum = exafs_spectrum + xas_spectrum.generate_from_contributions(logger) + if xas_spectrum.value is not None: + assert (xas_spectrum.value == xas_values).all() else: - assert xas_spectra.value == xas_values + assert xas_spectrum.value == xas_values diff --git a/tests/test_utils.py b/tests/test_utils.py index bdbe6e1c..5a498952 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -16,13 +16,20 @@ # limitations under the License. # -from . import logger +import pytest +import numpy as np + +from nomad.units import ureg from nomad_simulations.utils import ( get_sibling_section, is_not_representative, + get_variables, ) from nomad_simulations.model_system import ModelSystem, AtomicCell, Symmetry +from nomad_simulations.variables import Temperature, Energy2 as Energy + +from . import logger def test_get_sibling_section(): @@ -59,3 +66,28 @@ def test_is_not_representative(): # ! Missing test for RusselSandersState (but this class will probably be deprecated) + + +@pytest.mark.parametrize( + 'variables, result, result_length', + [ + (None, [], 0), + ([], [], 0), + ([Temperature()], [], 0), + ([Temperature(), Energy(n_points=4)], [Energy(n_points=4)], 1), + ( + [Temperature(), Energy(n_points=2), Energy(n_points=10)], + [Energy(n_points=2), Energy(n_points=10)], + 2, + ), + # TODO add testing when we have variables which inherit from another variable + ], +) +def test_get_variables(variables: list, result: list, result_length: int): + """ + Test the `get_variables` utility function + """ + energies = get_variables(variables, Energy) + assert len(energies) == result_length + for i, energy in enumerate(energies): # asserting energies == result does not work + assert energy.n_points == result[i].n_points