Add DielectricFunction properties (#72)

* Fix variables module and added Frequency

Commented out Variables.points.shape to fix when points are refs to a specific NumericalSettings section (like KMesh(Variables) or KLinePath(Variables))

Deleted refs to numerical settings under KMesh(Variables) and KLinePath(Variables)

* Changed to plural in the list of properties under Outputs

Added AbsorptionSpectrum to the list

* Added Permittivity property

Added extract_absorption_spectra from the imaginary part of the permittivity

Added AbsorptionSpectrum.axis to differentiate the principal axis to which is defined from

* Added get_variables utils function

Changed code to use this utils function

* Added testing for Permittivity and get_variables

Fixed testing

src/nomad_simulations/numerical_settings.py

@@ -336,7 +336,7 @@ class KMesh(Mesh):
         description="""
         Dictionary containing the high-symmetry point labels and their values in units of `reciprocal_lattice_vectors`.
         E.g., in a cubic lattice:
-            high_symmetry_points ={
+            high_symmetry_points = {
                 'Gamma': [0, 0, 0],
                 'X': [0.5, 0, 0],
                 'Y': [0, 0.5, 0],

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)
 
     # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
     # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

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

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

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)

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,
+)