Skip to content

Commit

Permalink
- Rewrite band_structure.py
Browse files Browse the repository at this point in the history 
- Turn `variables.py/KMesh` into a `Dataset`
  • Loading branch information
Nathan Daelman committed Nov 18, 2024
1 parent 8d1bf01 commit c318b7a
Show file tree
Hide file tree
Showing 2 changed files with 74 additions and 242 deletions.
281 changes: 62 additions & 219 deletions src/nomad_simulations/schema_packages/properties/band_structure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,84 +2,95 @@

import numpy as np
import pint
from nomad.config import config
from nomad.metainfo import Quantity, SubSection
from nomad.metainfo import MEnum, Quantity, SubSection
from nomad.metainfo.dataset import MDataset, Dataset
from nomad.datamodel.metainfo.physical_properties import PhysicalProperty
from ..variables import SpinChannel, KMesh

if TYPE_CHECKING:
from nomad.datamodel.datamodel import EntryArchive
from nomad.metainfo import Context, Section
from structlog.stdlib import BoundLogger

from nomad_simulations.schema_packages.atoms_state import AtomsState, OrbitalsState
from nomad_simulations.schema_packages.numerical_settings import KSpace
from nomad_simulations.schema_packages.physical_property import (
PhysicalProperty,
validate_quantity_wrt_value,
)
from nomad_simulations.schema_packages.properties.band_gap import ElectronicBandGap
from nomad_simulations.schema_packages.properties.fermi_surface import FermiSurface
from nomad_simulations.schema_packages.utils import get_sibling_section

from nomad.config import config
configuration = config.get_plugin_entry_point(
'nomad_simulations.schema_packages:nomad_simulations_plugin'
)

class ElectronicEigenvalues(MDataset):
m_def = PhysicalProperty(
type=np.float64,
unit='joule',
iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicEigenvalues',
description="""A base section used to define basic quantities for the `ElectronicEigenvalues` and `ElectronicBandStructure` properties.""",
default_variables=[SpinChannel, KMesh],
)

class BaseElectronicEigenvalues(PhysicalProperty):
"""
A base section used to define basic quantities for the `ElectronicEigenvalues` and `ElectronicBandStructure` properties.
"""

iri = ''
# ? Should we add functionalities to handle min/max of the `value` in some specific cases, .e.g, bands around the Fermi level,
# ? core bands separated by gaps, and equivalently, higher-energy valence bands separated by gaps?

n_bands = Quantity(
type=np.int32,
reciprocal_cell = Quantity(
type=KSpace.reciprocal_lattice_vectors,
description="""
Number of bands / eigenvalues.
Reference to the reciprocal lattice vectors stored under `KSpace`.
""",
) # !


class Occupancy(MDataset):
m_def = PhysicalProperty(
type=np.float64, # actually values between 0 - 2
iri = 'http://fairmat-nfdi.eu/taxonomy/Occupancy',
description="""
Electrons occupancy of an atom per orbital and spin. This is a number defined between 0 and 1 for
spin-polarized systems, and between 0 and 2 for non-spin-polarized systems. This property is
important when studying if an orbital or spin channel are fully occupied, at half-filling, or
fully emptied, which have an effect on the electron-electron interaction effects.
""",
default_variables=[SpinChannel, KMesh],
)

value = Quantity(
type=np.float64,
unit='joule',
atoms_state_ref = Quantity(
type=AtomsState,
description="""
Value of the electronic eigenvalues.
Reference to the `AtomsState` section in which the occupancy is calculated.
""",
)

def __init__(
self, m_def: 'Section' = None, m_context: 'Context' = None, **kwargs
) -> None:
super().__init__(m_def, m_context, **kwargs)
# ! `n_bands` need to be set up during initialization of the class
self.rank = [int(kwargs.get('n_bands'))]

def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
super().normalize(archive, logger)

orbitals_state_ref = Quantity(
type=OrbitalsState,
description="""
Reference to the `OrbitalsState` section in which the occupancy is calculated.
""",
)

class ElectronicEigenvalues(BaseElectronicEigenvalues):
""" """

iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicEigenvalues'
class ElectronicEigenstates(MDataset): # ? rename
m_def = Dataset(
type=bool,
default_variables=[ElectronicEigenvalues, Occupancy],
)

spin_channel = Quantity(
n_bands = Quantity(
type=np.int32,
description="""
Spin channel of the corresponding electronic eigenvalues. It can take values of 0 or 1.
Number of bands / eigenvalues.
""",
)

occupation = Quantity(
type=np.float64,
shape=['*', 'n_bands'],
kind = Quantity(
type=MEnum('KS', 'KSxc', 'SigX', 'SigC', 'Zk'),
description="""
Occupation of the electronic eigenvalues. This is a number depending whether the `spin_channel` has been set or not.
If `spin_channel` is set, then this number is between 0 and 1, where 0 means that the state is unoccupied and 1 means
that the state is fully occupied; if `spin_channel` is not set, then this number is between 0 and 2. The shape of
this quantity is defined as `[K.n_points, K.dimensionality, n_bands]`, where `K` is a `variable` which can
be `KMesh` or `KLinePath`, depending whether the simulation mapped the whole Brillouin zone or just a specific
path.
Contributions to the electronic eigenvalues. Example, in the case of a DFT+GW calculation, the GW eigenvalues
are stored under `value`, and each contribution is identified by `label`:
- `'KS'`: Kohn-Sham contribution. This is also stored in the DFT entry under `ElectronicEigenvalues.value`.
- `'KSxc'`: Diagonal matrix elements of the expectation value of the Kohn-Sahm exchange-correlation potential.
- `'SigX'`: Diagonal matrix elements of the exchange self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`.
- `'SigC'`: Diagonal matrix elements of the correlation self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`.
- `'Zk'`: Quasiparticle renormalization factors contribution. This is also stored in the GW entry under `QuasiparticleWeights.value`.
""",
)

Expand All @@ -101,45 +112,10 @@ class ElectronicEigenvalues(BaseElectronicEigenvalues):
""",
)

# ? Should we add functionalities to handle min/max of the `value` in some specific cases, .e.g, bands around the Fermi level,
# ? core bands separated by gaps, and equivalently, higher-energy valence bands separated by gaps?

value_contributions = SubSection(
sub_section=BaseElectronicEigenvalues.m_def,
repeats=True,
description="""
Contributions to the electronic eigenvalues. Example, in the case of a DFT+GW calculation, the GW eigenvalues
are stored under `value`, and each contribution is identified by `label`:
- `'KS'`: Kohn-Sham contribution. This is also stored in the DFT entry under `ElectronicEigenvalues.value`.
- `'KSxc'`: Diagonal matrix elements of the expectation value of the Kohn-Sahm exchange-correlation potential.
- `'SigX'`: Diagonal matrix elements of the exchange self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`.
- `'SigC'`: Diagonal matrix elements of the correlation self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`.
- `'Zk'`: Quasiparticle renormalization factors contribution. This is also stored in the GW entry under `QuasiparticleWeights.value`.
""",
)

reciprocal_cell = Quantity(
type=KSpace.reciprocal_lattice_vectors,
description="""
Reference to the reciprocal lattice vectors stored under `KSpace`.
""",
)

def __init__(
self, m_def: 'Section' = None, m_context: 'Context' = None, **kwargs
) -> None:
super().__init__(m_def, m_context, **kwargs)
self.name = self.m_def.name

@validate_quantity_wrt_value(name='occupation')
def order_eigenvalues(self) -> Union[bool, tuple[pint.Quantity, np.ndarray]]:
def order_eigenvalues(self) -> Union[bool, tuple['pint.Quantity', np.ndarray]]:
"""
Order the eigenvalues based on the `value` and `occupation`. The return `value` and
`occupation` are flattened.
Returns:
(Union[bool, tuple[pint.Quantity, np.ndarray]]): The flattened and sorted `value` and `occupation`. If validation
fails, then it returns `False`.
"""
total_shape = np.prod(self.value.shape)

Expand All @@ -160,15 +136,11 @@ def order_eigenvalues(self) -> Union[bool, tuple[pint.Quantity, np.ndarray]]:

def resolve_homo_lumo_eigenvalues(
self,
) -> tuple[Optional[pint.Quantity], Optional[pint.Quantity]]:
) -> tuple[Optional['pint.Quantity'], Optional['pint.Quantity']]:
"""
Resolve the `highest_occupied` and `lowest_unoccupied` eigenvalues by performing a binary search on the
flattened and sorted `value` and `occupation`. If these quantities already exist, overwrite them or return
them if it is not possible to resolve from `value` and `occupation`.
Returns:
(tuple[Optional[pint.Quantity], Optional[pint.Quantity]]): The `highest_occupied` and
`lowest_unoccupied` eigenvalues.
"""
# Sorting `value` and `occupation`
if not self.order_eigenvalues(): # validation fails
Expand Down Expand Up @@ -198,13 +170,10 @@ def resolve_homo_lumo_eigenvalues(

return homo, lumo

def extract_band_gap(self) -> Optional[ElectronicBandGap]:
def extract_band_gap(self) -> Optional['ElectronicBandGap']:
"""
Extract the electronic band gap from the `highest_occupied` and `lowest_unoccupied` eigenvalues.
If the difference of `highest_occupied` and `lowest_unoccupied` is negative, the band gap `value` is set to 0.0.
Returns:
(Optional[ElectronicBandGap]): The extracted electronic band gap section to be stored in `Outputs`.
"""
band_gap = None
homo, lumo = self.resolve_homo_lumo_eigenvalues()
Expand All @@ -217,59 +186,9 @@ def extract_band_gap(self) -> Optional[ElectronicBandGap]:
band_gap.value = lumo - homo
return band_gap

# TODO fix this method once `FermiSurface` property is implemented
def extract_fermi_surface(self, logger: 'BoundLogger') -> Optional[FermiSurface]:
"""
Extract the Fermi surface for metal systems and using the `FermiLevel.value`.
Args:
logger (BoundLogger): The logger to log messages.
Returns:
(Optional[FermiSurface]): The extracted Fermi surface section to be stored in `Outputs`.
"""
# Check if the system has a finite band gap
homo, lumo = self.resolve_homo_lumo_eigenvalues()
if (homo and lumo) and (lumo - homo).magnitude > 0:
return None

# Get the `fermi_level.value`
fermi_level = get_sibling_section(
section=self, sibling_section_name='fermi_level', logger=logger
)
if fermi_level is None:
logger.warning(
'Could not extract the `FermiSurface`, because `FermiLevel` is not stored.'
)
return None
fermi_level_value = fermi_level.value.magnitude

# Extract values close to the `fermi_level.value`
fermi_indices = np.logical_and(
self.value.magnitude
>= (fermi_level_value - configuration.fermi_surface_tolerance),
self.value.magnitude
<= (fermi_level_value + configuration.fermi_surface_tolerance),
)
fermi_values = self.value[fermi_indices]

# Store `FermiSurface` values
# ! This is wrong (!) the `value` should be the `KMesh.points`, not the `ElectronicEigenvalues.value`
fermi_surface = FermiSurface(
n_bands=self.n_bands,
is_derived=True,
physical_property_ref=self,
)
fermi_surface.variables = self.variables
fermi_surface.value = fermi_values
return fermi_surface

def resolve_reciprocal_cell(self) -> Optional[pint.Quantity]:
def resolve_reciprocal_cell(self) -> Optional['pint.Quantity']:
"""
Resolve the reciprocal cell from the `KSpace` numerical settings section.
Returns:
Optional[pint.Quantity]: _description_
"""
numerical_settings = self.m_xpath(
'm_parent.m_parent.model_method[-1].numerical_settings', dict=False
Expand Down Expand Up @@ -298,81 +217,5 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
if band_gap is not None:
self.m_parent.electronic_band_gaps.append(band_gap)

# TODO uncomment once `FermiSurface` property is implemented
# `FermiSurface` extraction
# fermi_surface = self.extract_fermi_surface(logger)
# if fermi_surface is not None:
# self.m_parent.fermi_surfaces.append(fermi_surface)

# Resolve `reciprocal_cell` from the `KSpace` numerical settings section
self.reciprocal_cell = self.resolve_reciprocal_cell()


class ElectronicBandStructure(ElectronicEigenvalues):
"""
Accessible energies by the charges (electrons and holes) in the reciprocal space.
"""

iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicBandStructure'

def __init__(
self, m_def: 'Section' = None, m_context: 'Context' = None, **kwargs
) -> None:
super().__init__(m_def, m_context, **kwargs)
self.name = self.m_def.name

def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
super().normalize(archive, logger)


class Occupancy(PhysicalProperty):
"""
Electrons occupancy of an atom per orbital and spin. This is a number defined between 0 and 1 for
spin-polarized systems, and between 0 and 2 for non-spin-polarized systems. This property is
important when studying if an orbital or spin channel are fully occupied, at half-filling, or
fully emptied, which have an effect on the electron-electron interaction effects.
"""

iri = 'http://fairmat-nfdi.eu/taxonomy/Occupancy'

atoms_state_ref = Quantity(
type=AtomsState,
description="""
Reference to the `AtomsState` section in which the occupancy is calculated.
""",
)

orbitals_state_ref = Quantity(
type=OrbitalsState,
description="""
Reference to the `OrbitalsState` section in which the occupancy is calculated.
""",
)

spin_channel = Quantity(
type=np.int32,
description="""
Spin channel of the corresponding electronic property. It can take values of 0 and 1.
""",
)

value = Quantity(
type=np.float64,
description="""
Value of the electronic occupancy in the atom defined by `atoms_state_ref` and the orbital
defined by `orbitals_state_ref`. the orbital. If `spin_channel` is set, then this number is
between 0 and 1, where 0 means that the state is unoccupied and 1 means that the state is
fully occupied; if `spin_channel` is not set, then this number is between 0 and 2.
""",
)

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

# TODO add extraction from `ElectronicEigenvalues.occupation`

def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
super().normalize(archive, logger)
self.reciprocal_cell = self.resolve_reciprocal_cell()
Loading

0 comments on commit c318b7a

Please sign in to comment.