diff --git a/properties/__init__.py b/properties/__init__.py new file mode 100644 index 00000000..c6b0a577 --- /dev/null +++ b/properties/__init__.py @@ -0,0 +1,32 @@ +# +# 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. + +from .energies import FermiLevel, ChemicalPotential +from .band_gap import ElectronicBandGap +from .spectral_profile import ( + SpectralProfile, + DOSProfile, + ElectronicDensityOfStates, + AbsorptionSpectrum, + XASSpectrum, +) +from .hopping_matrix import HoppingMatrix, CrystalFieldSplitting +from .permittivity import Permittivity +from .fermi_surface import FermiSurface +from .band_structure import ElectronicEigenvalues, ElectronicBandStructure + diff --git a/properties/band_gap.py b/properties/band_gap.py new file mode 100644 index 00000000..c6922a4f --- /dev/null +++ b/properties/band_gap.py @@ -0,0 +1,159 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import numpy as np +from structlog.stdlib import BoundLogger +import pint +from typing import Optional + +from nomad.metainfo import Quantity, MEnum, Section, Context + +from nomad_simulations.physical_property import PhysicalProperty + + +class ElectronicBandGap(PhysicalProperty): + """ + Energy difference between the highest occupied electronic state and the lowest unoccupied electronic state. + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicBandGap' + + type = Quantity( + type=MEnum('direct', 'indirect'), + description=""" + Type categorization of the electronic band gap. This quantity is directly related with `momentum_transfer` as by + definition, the electronic band gap is `'direct'` for zero momentum transfer (or if `momentum_transfer` is `None`) and `'indirect'` + for finite momentum transfer. + + Note: in the case of finite `variables`, this quantity refers to all of the `value` in the array. + """, + ) + + momentum_transfer = Quantity( + type=np.float64, + shape=[2, 3], + description=""" + If the electronic band gap is `'indirect'`, the reciprocal momentum transfer for which the band gap is defined + in units of the `reciprocal_lattice_vectors`. The initial and final momentum 3D vectors are given in the first + and second element. Example, the momentum transfer in bulk Si2 happens between the Γ and the (approximately) + X points in the Brillouin zone; thus: + `momentum_transfer = [[0, 0, 0], [0.5, 0.5, 0]]`. + + Note: this quantity only refers to scalar `value`, not to arrays of `value`. + """, + ) + + spin_channel = Quantity( + type=np.int32, + description=""" + Spin channel of the corresponding electronic band gap. It can take values of 0 or 1. + """, + ) + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the electronic band gap. This value has to be positive, otherwise it will + prop an error and be set to None by the `normalize()` function. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.name = self.m_def.name + self.rank = [] + + def validate_values(self, logger: BoundLogger) -> Optional[pint.Quantity]: + """ + Validate the electronic band gap `value` by checking if they are negative and sets them to None if they are. + + Args: + logger (BoundLogger): The logger to log messages. + """ + value = self.value.magnitude + if not isinstance(self.value.magnitude, np.ndarray): # for scalars + value = np.array( + [value] + ) # ! check this when talking with Lauri and Theodore + + # Set the value to 0 when it is negative + if (value < 0).any(): + logger.error('The electronic band gap cannot be defined negative.') + return None + + if not isinstance(self.value.magnitude, np.ndarray): # for scalars + value = value[0] + return value * self.value.u + + def resolve_type(self, logger: BoundLogger) -> Optional[str]: + """ + Resolves the `type` of the electronic band gap based on the stored `momentum_transfer` values. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[str]): The resolved `type` of the electronic band gap. + """ + mtr = self.momentum_transfer if self.momentum_transfer is not None else [] + + # Check if the `momentum_transfer` is [], and return the type and a warning in the log for `indirect` band gaps + if len(mtr) == 0: + if self.type == 'indirect': + logger.warning( + 'The `momentum_transfer` is not stored for an `indirect` band gap.' + ) + return self.type + + # Check if the `momentum_transfer` has at least two elements, and return None if it does not + if len(mtr) == 1: + logger.warning( + 'The `momentum_transfer` should have at least two elements so that the difference can be calculated and the type of electronic band gap can be resolved.' + ) + return None + + # Resolve `type` from the difference between the initial and final momentum transfer + momentum_difference = np.diff(mtr, axis=0) + if (np.isclose(momentum_difference, np.zeros(3))).all(): + return 'direct' + else: + return 'indirect' + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Checks if the `value` is negative and sets it to None if it is. + self.value = self.validate_values(logger) + if self.value is None: + # ? What about deleting the class if `value` is None? + logger.error('The `value` of the electronic band gap is not stored.') + return + + # Resolve the `type` of the electronic band gap from `momentum_transfer`, ONLY for scalar `value` + if isinstance(self.value.magnitude, np.ndarray): + logger.info( + 'We do not support `type` which describe individual elements in an array `value`.' + ) + else: + self.type = self.resolve_type(logger) + diff --git a/properties/band_structure.py b/properties/band_structure.py new file mode 100644 index 00000000..239503be --- /dev/null +++ b/properties/band_structure.py @@ -0,0 +1,328 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import numpy as np +from structlog.stdlib import BoundLogger +from typing import Optional, Tuple, Union +import pint + +from nomad.metainfo import Quantity, Section, Context, SubSection + +from nomad_simulations.numerical_settings import KSpace +from nomad_simulations.physical_property import ( + PhysicalProperty, + validate_quantity_wrt_value, +) +from nomad_simulations.properties import ElectronicBandGap, FermiSurface +from nomad_simulations.utils import get_sibling_section + + +class BaseElectronicEigenvalues(PhysicalProperty): + """ + A base section used to define basic quantities for the `ElectronicEigenvalues` and `ElectronicBandStructure` properties. + """ + + iri = '' + + n_bands = Quantity( + type=np.int32, + description=""" + Number of bands / eigenvalues. + """, + ) + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the electronic eigenvalues. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # ! `n_bands` need to be set up during initialization of the class + self.rank = [int(kwargs.get('n_bands'))] + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class ElectronicEigenvalues(BaseElectronicEigenvalues): + """ """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicEigenvalues' + + spin_channel = Quantity( + type=np.int32, + description=""" + Spin channel of the corresponding electronic eigenvalues. It can take values of 0 or 1. + """, + ) + + occupation = Quantity( + type=np.float64, + shape=['*', 'n_bands'], + description=""" + Occupation of the electronic eigenvalues. This is a number depending whether the `spin_channel` has been set or not. + If `spin_channel` is set, then this number is between 0 and 2, where 0 means that the state is unoccupied and 2 means + that the state is fully occupied; if `spin_channel` is not set, then this number is between 0 and 1. The shape of + this quantity is defined as `[K.n_points, K.dimensionality, n_bands]`, where `K` is a `variable` which can + be `KMesh` or `KLinePath`, depending whether the simulation mapped the whole Brillouin zone or just a specific + path. + """, + ) + + highest_occupied = Quantity( + type=np.float64, + unit='joule', + description=""" + Highest occupied electronic eigenvalue. Together with `lowest_unoccupied`, it defines the + electronic band gap. + """, + ) + + lowest_unoccupied = Quantity( + type=np.float64, + unit='joule', + description=""" + Lowest unoccupied electronic eigenvalue. Together with `highest_occupied`, it defines the + electronic band gap. + """, + ) + + # ? Should we add functionalities to handle min/max of the `value` in some specific cases, .e.g, bands around the Fermi level, + # ? core bands separated by gaps, and equivalently, higher-energy valence bands separated by gaps? + + value_contributions = SubSection( + sub_section=BaseElectronicEigenvalues.m_def, + repeats=True, + description=""" + Contributions to the electronic eigenvalues. Example, in the case of a DFT+GW calculation, the GW eigenvalues + are stored under `value`, and each contribution is identified by `label`: + - `'KS'`: Kohn-Sham contribution. This is also stored in the DFT entry under `ElectronicEigenvalues.value`. + - `'KSxc'`: Diagonal matrix elements of the expectation value of the Kohn-Sahm exchange-correlation potential. + - `'SigX'`: Diagonal matrix elements of the exchange self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`. + - `'SigC'`: Diagonal matrix elements of the correlation self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`. + - `'Zk'`: Quasiparticle renormalization factors contribution. This is also stored in the GW entry under `QuasiparticleWeights.value`. + """, + ) + + reciprocal_cell = Quantity( + type=KSpace.reciprocal_lattice_vectors, + description=""" + Reference to the reciprocal lattice vectors stored under `KSpace`. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.name = self.m_def.name + + @validate_quantity_wrt_value(name='occupation') + def order_eigenvalues(self) -> Union[bool, Tuple[pint.Quantity, np.ndarray]]: + """ + Order the eigenvalues based on the `value` and `occupation`. The return `value` and + `occupation` are flattened. + + Returns: + (Union[bool, Tuple[pint.Quantity, np.ndarray]]): The flattened and sorted `value` and `occupation`. If validation + fails, then it returns `False`. + """ + total_shape = np.prod(self.value.shape) + + # Order the indices in the flattened list of `value` + flattened_value = self.value.reshape(total_shape) + flattened_occupation = self.occupation.reshape(total_shape) + sorted_indices = np.argsort(flattened_value, axis=0) + + sorted_value = ( + np.take_along_axis(flattened_value.magnitude, sorted_indices, axis=0) + * flattened_value.u + ) + sorted_occupation = np.take_along_axis( + flattened_occupation, sorted_indices, axis=0 + ) + self.m_cache['sorted_eigenvalues'] = True + return sorted_value, sorted_occupation + + def resolve_homo_lumo_eigenvalues( + self, + ) -> Tuple[Optional[pint.Quantity], Optional[pint.Quantity]]: + """ + Resolve the `highest_occupied` and `lowest_unoccupied` eigenvalues by performing a binary search on the + flattened and sorted `value` and `occupation`. If these quantities already exist, overwrite them or return + them if it is not possible to resolve from `value` and `occupation`. + + Returns: + (Tuple[Optional[pint.Quantity], Optional[pint.Quantity]]): The `highest_occupied` and + `lowest_unoccupied` eigenvalues. + """ + # Sorting `value` and `occupation` + if not self.order_eigenvalues(): # validation fails + if self.highest_occupied is not None and self.lowest_unoccupied is not None: + return self.highest_occupied, self.lowest_unoccupied + return None, None + sorted_value, sorted_occupation = self.order_eigenvalues() + sorted_value_unit = sorted_value.u + sorted_value = sorted_value.magnitude + + # Binary search ot find the transition point between `occupation = 2` and `occupation = 0` + tolerance = 1e-3 # TODO add tolerance from config fields + homo = self.highest_occupied + lumo = self.lowest_unoccupied + mid = np.searchsorted(sorted_occupation <= tolerance, True) - 1 + if mid >= 0 and mid < len(sorted_occupation) - 1: + if sorted_occupation[mid] > 0 and ( + sorted_occupation[mid + 1] >= -tolerance + and sorted_occupation[mid + 1] <= tolerance + ): + homo = sorted_value[mid] * sorted_value_unit + lumo = sorted_value[mid + 1] * sorted_value_unit + + return homo, lumo + + def extract_band_gap(self) -> Optional[ElectronicBandGap]: + """ + Extract the electronic band gap from the `highest_occupied` and `lowest_unoccupied` eigenvalues. + If the difference of `highest_occupied` and `lowest_unoccupied` is negative, the band gap `value` is set to 0.0. + + Returns: + (Optional[ElectronicBandGap]): The extracted electronic band gap section to be stored in `Outputs`. + """ + band_gap = None + homo, lumo = self.resolve_homo_lumo_eigenvalues() + if homo and lumo: + band_gap = ElectronicBandGap(is_derived=True, physical_property_ref=self) + + if (lumo - homo).magnitude < 0: + band_gap.value = 0.0 + else: + band_gap.value = lumo - homo + return band_gap + + # TODO fix this method once `FermiSurface` property is implemented + def extract_fermi_surface(self, logger: BoundLogger) -> Optional[FermiSurface]: + """ + Extract the Fermi surface for metal systems and using the `FermiLevel.value`. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[FermiSurface]): The extracted Fermi surface section to be stored in `Outputs`. + """ + # Check if the system has a finite band gap + homo, lumo = self.resolve_homo_lumo_eigenvalues() + if (homo and lumo) and (lumo - homo).magnitude > 0: + return None + + # Get the `fermi_level.value` + fermi_level = get_sibling_section( + section=self, sibling_section_name='fermi_level', logger=logger + ) + if fermi_level is None: + logger.warning( + 'Could not extract the `FermiSurface`, because `FermiLevel` is not stored.' + ) + return None + fermi_level_value = fermi_level.value.magnitude + + # Extract values close to the `fermi_level.value` + tolerance = 1e-8 # TODO change this for a config field + fermi_indices = np.logical_and( + self.value.magnitude >= (fermi_level_value - tolerance), + self.value.magnitude <= (fermi_level_value + tolerance), + ) + fermi_values = self.value[fermi_indices] + + # Store `FermiSurface` values + # ! This is wrong (!) the `value` should be the `KMesh.points`, not the `ElectronicEigenvalues.value` + fermi_surface = FermiSurface( + n_bands=self.n_bands, + is_derived=True, + physical_property_ref=self, + ) + fermi_surface.variables = self.variables + fermi_surface.value = fermi_values + return fermi_surface + + def resolve_reciprocal_cell(self) -> Optional[pint.Quantity]: + """ + Resolve the reciprocal cell from the `KSpace` numerical settings section. + + Returns: + Optional[pint.Quantity]: _description_ + """ + numerical_settings = self.m_xpath( + 'm_parent.m_parent.model_method[-1].numerical_settings', dict=False + ) + if numerical_settings is None: + return None + k_space = None + for setting in numerical_settings: + if isinstance(setting, KSpace): + k_space = setting + break + if k_space is None: + return None + return k_space + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Resolve `highest_occupied` and `lowest_unoccupied` eigenvalues + self.highest_occupied, self.lowest_unoccupied = ( + self.resolve_homo_lumo_eigenvalues() + ) + + # `ElectronicBandGap` extraction + band_gap = self.extract_band_gap() + if band_gap is not None: + self.m_parent.electronic_band_gaps.append(band_gap) + + # TODO uncomment once `FermiSurface` property is implemented + # `FermiSurface` extraction + # fermi_surface = self.extract_fermi_surface(logger) + # if fermi_surface is not None: + # self.m_parent.fermi_surfaces.append(fermi_surface) + + # Resolve `reciprocal_cell` from the `KSpace` numerical settings section + self.reciprocal_cell = self.resolve_reciprocal_cell() + + +class ElectronicBandStructure(ElectronicEigenvalues): + """ + Accessible energies by the charges (electrons and holes) in the reciprocal space. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicBandStructure' + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + diff --git a/properties/energies.py b/properties/energies.py new file mode 100644 index 00000000..8ef7c7fb --- /dev/null +++ b/properties/energies.py @@ -0,0 +1,801 @@ +# +# 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. +# +# +# 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, SubSection, MEnum +from .physical_property import PhysicalProperty + +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import numpy as np + +from nomad.metainfo import Quantity, Section, Context + +from nomad_simulations.physical_property import PhysicalProperty + + +class FermiLevel(PhysicalProperty): + """ + Energy required to add or extract a charge from a material at zero temperature. It can be also defined as the chemical potential at zero temperature. + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + iri = 'http://fairmat-nfdi.eu/taxonomy/FermiLevel' + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the Fermi level. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.rank = [] + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class ChemicalPotential(PhysicalProperty): + """ + Free energy cost of adding or extracting a particle from a thermodynamic system. + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + iri = 'http://fairmat-nfdi.eu/taxonomy/ChemicalPotential' + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the chemical potential. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.rank = [] + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class EnergyContributions(PhysicalProperty): + """ + Base section for all contribution groups to the total energy. + """ + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +#################################################### +# List of classical energy contribuions +#################################################### + +class PotentialEnergy(PhysicalProperty): + """ + Section containing the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class KineticEnergy(PhysicalProperty): + """ + Section containing the kinetic energy of a (sub)system. + """ + + type = Quantity( + type=MEnum('classical', 'quantum'), + description=""" + Refers to the method used for calculating the kinetic energy. + + Allowed values are: + + | Thermostat Name | Description | + + | ---------------------- | ----------------------------------------- | + + | `"classical"` | The kinetic energy is calculated directly from + the velocities of particles as KE = /sum_i 1/2 m_i v_i^2, where the sum runs over the number of particles in the system, + m_i is the mass of particle i, and v_i is the velocity of particle i. | + + | `"quantum"` | ... | + """, + ) + # ? Just an idea, not sure if this is necessary since we will be referencing the method, I guess maybe it's overkill in general + # ? If we do this, maybe we should have a general Energy(PP) class which is inherited from in case there are multiple interpretations/calculation methods + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the kinetic energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class VDWEnergy(PhysicalProperty): + """ + Section containing the van der Waals contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the van der Waals contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectrostaticEnergy(PhysicalProperty): + """ + Section containing all electrostatic contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of all electrostatic contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectrostaticShortRangeEnergy(PhysicalProperty): + """ + Section containing short-range electrostatic contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the short-range electrostatic contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectrostaticLongRangeEnergy(PhysicalProperty): + """ + Section containing long-range electrostatic contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the long-range electrostatic contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class BondedEnergy(PhysicalProperty): + """ + Section containing all bonded (i.e., intramolecular) contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of all bonded (i.e., intramolecular) contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class BondEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from bond interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from bond interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class AngleEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from angle interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from angle interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class DihedralEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from dihedral interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from dihedral interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ImproperDihedralEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from improper dihedral interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from improper dihedral interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ExternalEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from external interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from external interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class EnergyContributionsClassical(EnergyContributions): + """ + Section containing contributions to the potential energy from a classical force field. + """ + + potential = SubSection(sub_section=PotentialEnergy.m_def, repeats=False) + + kinetic = SubSection(sub_section=KineticEnergy.m_def, repeats=False) + + vdw = SubSection(sub_section=VDWEnergy.m_def, repeats=False) + + electrostatic = SubSection(sub_section=ElectrostaticEnergy.m_def, repeats=False) + + electrostatic_short_range = SubSection(sub_section=ElectrostaticShortRangeEnergy.m_def, repeats=False) + + electrostatic_long_range = SubSection(sub_section=ElectrostaticLongRangeEnergy.m_def, repeats=False) + + bonded = SubSection(sub_section=BondedEnergy.m_def, repeats=False) + + bond = SubSection(sub_section=BondEnergy.m_def, repeats=False) + + angle = SubSection(sub_section=AngleEnergy.m_def, repeats=False) + + dihedral = SubSection(sub_section=DihedralEnergy.m_def, repeats=False) + + improper_dihedral = SubSection(sub_section=ImproperDihedralEnergy.m_def, repeats=False) + + external = SubSection(sub_section=ExternalEnergy.m_def, repeats=False) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + # Set the name of the section + self.name = self.m_def.name + +###################################### +# List of quantum energy contributions +###################################### + +class ElectronicEnergy(PhysicalProperty): + """ + Section containing the electronic energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the electronic energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectronicKineticEnergy(PhysicalProperty): + """ + Section containing the electronic kinetic energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the electronic kinetic energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class XCEnergy(PhysicalProperty): + """ + Section containing the exchange-correlation (XC) energy of a (sub)system, + calculated using the functional stored in XC_functional. + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the exchange-correlation energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class XCPotentialEnergy(PhysicalProperty): + """ + Section containing the potential energy contribution to the exchange-correlation (XC) energy, + i.e., the integral of the first order derivative of the functional + stored in XC_functional (integral of v_xc*electron_density), i.e., the component + of XC that is in the sum of the eigenvalues. Value associated with the + configuration, should be the most converged value. | + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the potential energy contribution of the exchange-correlation energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +# ? XCCorrelationEnergy? +class CorrelationEnergy(PhysicalProperty): + """ + Section containing the correlation energy of a (sub)system, + calculated using the method described in XC_functional. + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the correlation energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +# ? XCExchangeEnergy? +class ExchangeEnergy(PhysicalProperty): + """ + Section containing the exchange energy of a (sub)system, + calculated using the method described in XC_functional. + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the exchange energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class T0Energy(PhysicalProperty): + """ + Section containing the total energy of a (sub)system extrapolated to $T=0$. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the total energy extrapolated to + $T=0$, based on a free-electron gas argument. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ZeroPointEnergy(PhysicalProperty): + """ + Section containing the zero-point vibrational energy of a (sub)system, + calculated using the method described in zero_point_method. + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the zero-point energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectrostaticEnergy(PhysicalProperty): + """ + Section containing the electrostatic energy (nuclei + electrons) of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the electrostatic energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class NuclearRepulsionEnergy(PhysicalProperty): + """ + Section containing the nuclear-nuclear repulsion energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the nuclear-nuclear repulsion energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class EnergyTotalContributionsQuantum(EnergyContributions): + """ + Section containing contributions to the potential energy from a DFT calculation. + """ +... + + electronic = SubSection(sub_section=ElectronicEnergy.m_def, repeats=False) + + electronic_kinetic = SubSection(sub_section=ElectronicKineticEnergy.m_def, repeats=False) + + xc = SubSection(sub_section=XCEnergy.m_def, repeats=False) + + xc_potential = SubSection(sub_section=XCPotentialEnergy.m_def, repeats=False) + + correlation = SubSection(sub_section=CorrelationEnergy.m_def, repeats=False) + + exchange = SubSection(sub_section=ExchangeEnergy.m_def, repeats=False) + + t0 = SubSection(sub_section=T0Energy.m_def, repeats=False) + + zero_point = SubSection(sub_section=ZeroPointEnergy.m_def, repeats=False) + + electrostatic = SubSection(sub_section=ElectrostaticEnergy.m_def, repeats=False) + + nuclear_repulsion = SubSection(sub_section=NuclearRepulsionEnergy.m_def, repeats=False) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +########################## +# Other / General energies +########################## + +class EnergyTotal(PhysicalProperty): + """ + Section containing the total energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the total energy. + """, + ) + # ? Do we need these descriptions under value? It ends up simply duplicating the section info to some extent. + + # ? I think we want to allow some flexible to the contributions here, like so: ?? + contributions = SubSection(sub_section=EnergyContributions.m_def, repeats=True) + # classical_contributions = SubSection(sub_section=EnergyContributionsClassical.m_def, repeats=False) + # quantum_contributions = SubSection(sub_section=EnergyTotalContributionsQuantum.m_def, repeats=False) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # madelung = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the Madelung energy. + # """, + # ) + + # # TODO I suggest ewald is moved to "long range" under electrostatic->energyentry, unless there is some other usage I am misunderstanding + # ewald = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the Ewald energy. + # """, + # ) + + # free = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the free energy (nuclei + electrons) + # (whose minimum gives the smeared occupation density calculated with + # smearing_kind). + # """, + # ) + + # sum_eigenvalues = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the sum of the eigenvalues of the + # Hamiltonian matrix. + # """, + # ) + + # van_der_waals = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the Van der Waals energy. A multiple + # occurence is expected when more than one van der Waals methods are defined. The + # van der Waals kind should be specified in Energy.kind + # """, + # ) + + # hartree_fock_x_scaled = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Scaled exact-exchange energy that depends on the mixing parameter of the + # functional. For example in hybrid functionals, the exchange energy is given as a + # linear combination of exact-energy and exchange energy of an approximate DFT + # functional; the exact exchange energy multiplied by the mixing coefficient of the + # hybrid functional would be stored in this metadata. Defined consistently with + # XC_method. + # """, + # ) + + # ? This is technically NOT redundant with total energy, but I fear that people will use them interchangeably, so we need to be clear about the definitions in any case + # internal = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Value of the internal energy. + # """, + # ) + + # double_counting = SubSection( + # sub_section=EnergyEntry.m_def, + # categories=[FastAccess], + # description=""" + # Double counting correction when performing Hubbard model calculations. + # """, + # ) + + # # TODO remove this should be be entropy.correction + # correction_entropy = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Entropy correction to the potential energy to compensate for the change in + # occupation so that forces at finite T do not need to keep the change of occupation + # in account. Defined consistently with XC_method. + # """, + # ) + + # # TODO remove this should be in electrostatic.correction + # correction_hartree = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Correction to the density-density electrostatic energy in the sum of eigenvalues + # (that uses the mixed density on one side), and the fully consistent density- + # density electrostatic energy. Defined consistently with XC_method. + # """, + # ) + + # # TODO remove this should be in xc.correction + # correction_xc = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Correction to energy_XC. + # """, + # ) + + # # ? Is it ok to store this in outputs and not in workflow? I guess we can ensure in normalization that this is a WorkflowOutput, etc...? + # change = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Stores the change of total energy with respect to the previous step. + # """, + # categories=[ErrorEstimateContribution, EnergyValue], + # ) + + # fermi = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Fermi energy (separates occupied from unoccupied single-particle states) + # """, + # categories=[EnergyTypeReference, EnergyValue], + # ) + + # highest_occupied = Quantity( + # type=np.dtype(np.float64), + # unit='joule', + # shape=[], + # description=""" + # The highest occupied energy. + # """, + # ) + + # lowest_unoccupied = Quantity( + # type=np.dtype(np.float64), + # unit='joule', + # shape=[], + # description=""" + # The lowest unoccupied energy. + # """, + # ) + + + # # TODO this should be removed and replaced by correction in EnergyEntry + # current = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the energy calculated with + # calculation_method_current. energy_current is equal to energy_total for + # non-perturbative methods. For perturbative methods, energy_current is equal to the + # correction: energy_total minus energy_total of the calculation_to_calculation_ref + # with calculation_to_calculation_kind = starting_point + # """, + # ) + + + + +# # ? Do we want to allow this? +# class EnergyCustom(PhysicalProperty): +# """ +# Section containing the total energy of a (sub)system. +# """ + +# name = Quantity( +# type=str, +# description=""" +# The name of this custom energy. +# """, +# ) + +# value = Quantity( +# type=np.float64, +# unit='joule', +# description=""" +# The value of this custom energy type. +# """, +# ) + +# def normalize(self, archive, logger) -> None: +# super().normalize(archive, logger) + diff --git a/properties/fermi_surface.py b/properties/fermi_surface.py new file mode 100644 index 00000000..4b473a9a --- /dev/null +++ b/properties/fermi_surface.py @@ -0,0 +1,53 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import numpy as np + +from nomad.metainfo import Quantity, Section, Context + +from nomad_simulations.physical_property import PhysicalProperty + + +# TODO This class is not implemented yet. @JosePizarro3 will work in another PR to implement it. +class FermiSurface(PhysicalProperty): + """ + Energy boundary in reciprocal space that separates the filled and empty electronic states in a metal. + It is related with the crossing points in reciprocal space by the chemical potential or, equivalently at + zero temperature, the Fermi level. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/FermiSurface' + + n_bands = Quantity( + type=np.int32, + description=""" + Number of bands / eigenvalues. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # ! `n_bands` need to be set up during initialization of the class + self.rank = [int(kwargs.get('n_bands'))] + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + diff --git a/properties/forces.py b/properties/forces.py new file mode 100644 index 00000000..a34f4508 --- /dev/null +++ b/properties/forces.py @@ -0,0 +1,182 @@ +# +# 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. +# +# +# 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 re +import numpy as np + +from nomad.units import ureg +from nomad.datamodel.data import ArchiveSection +from nomad.metainfo import Quantity, SubSection, SectionProxy, MEnum +from nomad.metainfo.metainfo import DirectQuantity, Dimension +from nomad.datamodel.metainfo.annotations import ELNAnnotation + +from .outputs import Outputs, SCFOutputs, WorkflowOutputs, TrajectoryOutputs +from .atoms_state import AtomsState +from .physical_property import PhysicalProperty + + +from nomad.datamodel.metainfo.common import ( + FastAccess, + PropertySection, + ProvenanceTracker, +) +from nomad.metainfo import ( + Category, + HDF5Reference, + MCategory, + MEnum, + MSection, + Package, + Quantity, + Reference, + Section, + SectionProxy, + SubSection, +) + +from .model_system import ModelSystem + + + +class ForceTotal(PhysicalProperty): + """ + Section containing the total force of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='newton', + description=""" + The value of the total force. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + + +class ForcesEntry(Atomic): + """ + Section describing a contribution to or type of atomic forces. + """ + + m_def = Section(validate=False) + + value = Quantity( + type=np.dtype(np.float64), + shape=['n_atoms', 3], + unit='newton', + description=""" + Value of the forces acting on the atoms. This is calculated as minus gradient of + the corresponding energy type or contribution **including** constraints, if + present. The derivatives with respect to displacements of nuclei are evaluated in + Cartesian coordinates. In addition, these are obtained by filtering out the + unitary transformations (center-of-mass translations and rigid rotations for + non-periodic systems, see value_raw for the unfiltered counterpart). + """, + ) + + value_raw = Quantity( + type=np.dtype(np.float64), + shape=['n_atoms', 3], + unit='newton', + description=""" + Value of the forces acting on the atoms **not including** such as fixed atoms, + distances, angles, dihedrals, etc.""", + ) + # ? This is VERY imprecise, is this used regularly? + +class Forces(MSection): + """ + Section containing all forces types and contributions. + """ + + m_def = Section(validate=False) + + total = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains the value and information regarding the total forces on the atoms + calculated as minus gradient of energy_total. + """, + ) + # ! We need to avoid giving the precise method of calculation without also providing context, this is not necessarily true in general! + + free = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains the value and information regarding the forces on the atoms + corresponding to the minus gradient of energy_free. The (electronic) energy_free + contains the information on the change in (fractional) occupation of the + electronic eigenstates, which are accounted for in the derivatives, yielding a + truly energy-conserved quantity. + """, + ) + + t0 = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains the value and information regarding the forces on the atoms + corresponding to the minus gradient of energy_T0. + """, + ) + + contributions = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains other forces contributions to the total atomic forces not already + defined. + """, + repeats=True, + ) + + types = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains other types of forces not already defined. + """, + repeats=True, + ) + + + + + + forces = SubSection(sub_section=Forces.m_def) + + + diff --git a/properties/hopping_matrix.py b/properties/hopping_matrix.py new file mode 100644 index 00000000..cebeea71 --- /dev/null +++ b/properties/hopping_matrix.py @@ -0,0 +1,105 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import numpy as np + +from nomad.metainfo import Quantity, Section, Context + +from nomad_simulations.physical_property import PhysicalProperty + + +class HoppingMatrix(PhysicalProperty): + """ + Transition probability between two atomic orbitals in a tight-binding model. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/HoppingMatrix' + + n_orbitals = Quantity( + type=np.int32, + description=""" + Number of orbitals in the tight-binding model. + """, + ) + + degeneracy_factors = Quantity( + type=np.int32, + shape=['*'], + description=""" + Degeneracy of each Wigner-Seitz point. + """, + ) + + value = Quantity( + type=np.complex128, + unit='joule', + description=""" + Value of the hopping matrix in joules. The elements are complex numbers defined for each Wigner-Seitz point and + each pair of orbitals; thus, `rank = [n_orbitals, n_orbitals]`. Note this contains also the onsite values, i.e., + it includes the Wigner-Seitz point (0, 0, 0), hence the `CrystalFieldSplitting` values. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # ! n_orbitals need to be set up during initialization of the class + self.rank = [int(kwargs.get('n_orbitals')), int(kwargs.get('n_orbitals'))] + self.name = self.m_def.name + + # TODO add normalization to extract DOS, band structure, etc, properties from `HoppingMatrix` + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class CrystalFieldSplitting(PhysicalProperty): + """ + Energy difference between the degenerated orbitals of an ion in a crystal field environment. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/CrystalFieldSplitting' + + n_orbitals = Quantity( + type=np.int32, + description=""" + Number of orbitals in the tight-binding model. + """, + ) + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the crystal field splittings in joules. This is the intra-orbital local contribution, i.e., the same orbital + at the same Wigner-Seitz point (0, 0, 0). + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # ! `n_orbitals` need to be set up during initialization of the class + self.rank = [int(kwargs.get('n_orbitals'))] + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + diff --git a/properties/permittivity.py b/properties/permittivity.py new file mode 100644 index 00000000..d97c5cd4 --- /dev/null +++ b/properties/permittivity.py @@ -0,0 +1,119 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import numpy as np +from structlog.stdlib import BoundLogger +from typing import Optional, List + +from nomad.metainfo import Quantity, Section, Context, MEnum + +from nomad_simulations.physical_property import PhysicalProperty +from nomad_simulations.utils import get_variables +from nomad_simulations.variables import Frequency, KMesh +from nomad_simulations.properties.spectral_profile import AbsorptionSpectrum + + +# TODO add `DielectricStrength` when we have examples and understand how to extract it from the `Permittivity` tensor. + + +class Permittivity(PhysicalProperty): + """ + Response of the material to polarize in the presence of an electric field. + + Alternative names: `DielectricFunction`. + """ + + iri = 'http://fairmat-nfdi.eu/taxonomy/Permittivity' + + type = Quantity( + type=MEnum('static', 'dynamic'), + description=""" + Type of permittivity which allows to identify if the permittivity depends on the frequency or not. + """, + ) + + value = Quantity( + type=np.complex128, + # unit='joule', # TODO check units (they have to match `SpectralProfile.value`) + description=""" + Value of the permittivity tensor. If the value does not depend on the scattering vector `q`, then we + can extract the optical absorption spectrum from the imaginary part of the permittivity tensor (this is also called + macroscopic dielectric function). + """, + ) + + # ? We need use cases to understand if we need to define contributions to the permittivity tensor. + # ? `ionic` and `electronic` contributions are common in the literature. + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.rank = [3, 3] + self.name = self.m_def.name + self._axes_map = ['xx', 'yy', 'zz'] + + def resolve_type(self) -> str: + frequencies = get_variables(self.variables, Frequency) + if len(frequencies) > 0: + return 'dynamic' + return 'static' + + def extract_absorption_spectra( + self, logger: BoundLogger + ) -> Optional[List[AbsorptionSpectrum]]: + """ + Extract the absorption spectrum from the imaginary part of the permittivity tensor. + """ + # If the `pemittivity` depends on the scattering vector `q`, then we cannot extract the absorption spectrum + q_mesh = get_variables(self.variables, KMesh) + if len(q_mesh) > 0: + logger.warning( + 'The `permittivity` depends on the scattering vector `q`, so that we cannot extract the absorption spectrum.' + ) + return None + # Extract the `Frequency` variable to extract the absorption spectrum + frequencies = get_variables(self.variables, Frequency) + if len(frequencies) == 0: + logger.warning( + 'The `permittivity` does not have a `Frequency` variable to extract the absorption spectrum.' + ) + return None + # Define the `absorption_spectra` for each principal direction along the diagonal of the `Permittivity.value` as the imaginary part + spectra = [] + for i in range(3): + val = self.value[:, i, i].imag + absorption_spectrum = AbsorptionSpectrum( + axis=self._axes_map[i], variables=frequencies + ) + absorption_spectrum.value = val + absorption_spectrum.physical_property_ref = self + spectra.append(absorption_spectrum) + return spectra + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Resolve the `type` of permittivity + self.type = self.resolve_type() + + # `AbsorptionSpectrum` extraction + absorption_spectra = self.extract_absorption_spectra(logger) + if absorption_spectra is not None: + self.m_parent.absorption_spectrum = absorption_spectra + diff --git a/properties/spectral_profile.py b/properties/spectral_profile.py new file mode 100644 index 00000000..f2c894f0 --- /dev/null +++ b/properties/spectral_profile.py @@ -0,0 +1,619 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import numpy as np +from structlog.stdlib import BoundLogger +from typing import Optional, List, Dict +import pint + +from nomad import config +from nomad.metainfo import Quantity, SubSection, Section, Context, MEnum + +from nomad_simulations.utils import get_sibling_section, get_variables +from nomad_simulations.physical_property import PhysicalProperty +from nomad_simulations.variables import Energy2 as Energy +from nomad_simulations.atoms_state import AtomsState, OrbitalsState +from nomad_simulations.properties.band_gap import ElectronicBandGap + + +class SpectralProfile(PhysicalProperty): + """ + A base section used to define the spectral profile. + """ + + value = Quantity( + type=np.float64, + description=""" + The value of the intensities of a spectral profile in arbitrary units. + """, + ) # TODO check units and normalization_factor of DOS and Spectras and see whether they can be merged + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.rank = [] + + def is_valid_spectral_profile(self) -> bool: + """ + Check if the spectral profile is valid, i.e., if all `value` are defined positive. + + Returns: + (bool): True if the spectral profile is valid, False otherwise. + """ + if (self.value < 0.0).any(): + return False + return True + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + if self.is_valid_spectral_profile() is False: + logger.error( + 'Invalid negative intensities found: could not validate spectral profile.' + ) + return + + +class DOSProfile(SpectralProfile): + """ + A base section used to define the `value` of the `ElectronicDensityOfState` property. This is useful when containing + contributions for `projected_dos` with the correct unit. + """ + + value = Quantity( + type=np.float64, + unit='1/joule', + description=""" + The value of the electronic DOS. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + + def resolve_pdos_name(self, logger: BoundLogger) -> Optional[str]: + """ + Resolve the `name` of the projected `DOSProfile` from the `entity_ref` section. This is resolved as: + - `'atom X'` with 'X' being the chemical symbol for `AtomsState` references. + - `'orbital Y X'` with 'X' being the chemical symbol and 'Y' the orbital label for `OrbitalsState` references. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[str]): The resolved `name` of the projected DOS profile. + """ + if self.entity_ref is None: + logger.warning( + 'The `entity_ref` is not set for the DOS profile. Could not resolve the `name`.' + ) + return None + + # Resolve the `name` from the `entity_ref` + name = None + if isinstance(self.entity_ref, AtomsState): + name = f'atom {self.entity_ref.chemical_symbol}' + elif isinstance(self.entity_ref, OrbitalsState): + name = f'orbital {self.entity_ref.l_quantum_symbol}{self.entity_ref.ml_quantum_symbol} {self.entity_ref.m_parent.chemical_symbol}' + return name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # We resolve + self.name = self.resolve_pdos_name(logger) + + +class ElectronicDensityOfStates(DOSProfile): + """ + Number of electronic states accessible for the charges per energy and per volume. + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicDensityOfStates' + + spin_channel = Quantity( + type=np.int32, + description=""" + Spin channel of the corresponding electronic DOS. It can take values of 0 or 1. + """, + ) + + # TODO clarify the role of `energies_origin` once `ElectronicEigenvalues` is implemented + energies_origin = Quantity( + type=np.float64, + unit='joule', + description=""" + Energy level denoting the origin along the energy axis, used for comparison and visualization. It is + defined as the `ElectronicEigenvalues.highest_occupied_energy`. + """, + ) + + normalization_factor = Quantity( + type=np.float64, + description=""" + Normalization factor for electronic DOS to get a cell-independent intensive DOS. The cell-independent + intensive DOS is as the integral from the lowest (most negative) energy to the Fermi level for a neutrally + charged system (i.e., the sum of `AtomsState.charge` is zero). + """, + ) + + # ? Do we want to store the integrated value here os as part of an nomad-analysis tool? Check `dos_integrator.py` module in dos normalizer repository + # value_integrated = Quantity( + # type=np.float64, + # description=""" + # The cumulative intensities integrated from from the lowest (most negative) energy to the Fermi level. + # """, + # ) + + projected_dos = SubSection( + sub_section=DOSProfile.m_def, + repeats=True, + description=""" + Projected DOS. It can be atom- (different elements in the unit cell) or orbital-projected. These can be calculated in a cascade as: + - If the total DOS is not present, we sum all atom-projected DOS to obtain it. + - If the atom-projected DOS is not present, we sum all orbital-projected DOS to obtain it. + Note: the cover given by summing up contributions is not perfect, and will depend on the projection functions used. + + In `projected_dos`, `name` and `entity_ref` must be set in order for normalization to work: + - The `entity_ref` is the `OrbitalsState` or `AtomsState` sections. + - The `name` of the projected DOS should be `'atom X'` or `'orbital Y X'`, with 'X' being the chemical symbol and 'Y' the orbital label. + These can be extracted from `entity_ref`. + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + self.name = self.m_def.name + + def resolve_energies_origin( + self, + energies: pint.Quantity, + fermi_level: Optional[pint.Quantity], + logger: BoundLogger, + ) -> Optional[pint.Quantity]: + """ + Resolve the origin of reference for the energies from the sibling `ElectronicEigenvalues` section and its + `highest_occupied` level, or if this does not exist, from the `fermi_level` value as extracted from the sibling property, `FermiLevel`. + + Args: + fermi_level (Optional[pint.Quantity]): The resolved Fermi level. + energies (pint.Quantity): The grid points of the `Energy` variable. + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[pint.Quantity]): The resolved origin of reference for the energies. + """ + # Check if the variables contain more than one variable (different than Energy) + # ? Is this correct or should be use the index of energies to extract the proper shape element in `self.value` being used for `dos_values`? + if len(self.variables) > 1: + logger.warning( + 'The ElectronicDensityOfStates section contains more than one variable. We cannot extract the energy reference.' + ) + return None + + # Extract the `ElectronicEigenvalues` section to get the `highest_occupied` and `lowest_unoccupied` energies + # TODO implement once `ElectronicEigenvalues` is in the schema + eigenvalues = get_sibling_section( + section=self, sibling_section_name='electronic_eigenvalues', logger=logger + ) # we consider `index_sibling` to be 0 + highest_occupied_energy = ( + eigenvalues.highest_occupied if eigenvalues is not None else None + ) + lowest_unoccupied_energy = ( + eigenvalues.lowest_unoccupied if eigenvalues is not None else None + ) + # and set defaults for `highest_occupied_energy` and `lowest_unoccupied_energy` in `m_cache` + if highest_occupied_energy is not None: + self.m_cache['highest_occupied_energy'] = highest_occupied_energy + if lowest_unoccupied_energy is not None: + self.m_cache['lowest_unoccupied_energy'] = lowest_unoccupied_energy + + # Set thresholds for the energies and values + energy_threshold = config.normalize.band_structure_energy_tolerance + value_threshold = 1e-8 # The DOS value that is considered to be zero + + # Check that the closest `energies` to the energy reference is not too far away. + # If it is very far away, normalization may be very inaccurate and we do not report it. + dos_values = self.value.magnitude + eref = highest_occupied_energy if fermi_level is None else fermi_level + fermi_idx = (np.abs(energies - eref)).argmin() + fermi_energy_closest = energies[fermi_idx] + distance = np.abs(fermi_energy_closest - eref) + single_peak_fermi = False + if distance.magnitude <= energy_threshold: + # See if there are zero values close below the energy reference. + idx = fermi_idx + idx_descend = fermi_idx + while True: + try: + value = dos_values[idx] + energy_distance = np.abs(eref - energies[idx]) + except IndexError: + break + if energy_distance.magnitude > energy_threshold: + break + if value <= value_threshold: + idx_descend = idx + break + idx -= 1 + + # See if there are zero values close above the fermi energy. + idx = fermi_idx + idx_ascend = fermi_idx + while True: + try: + value = dos_values[idx] + energy_distance = np.abs(eref - energies[idx]) + except IndexError: + break + if energy_distance.magnitude > energy_threshold: + break + if value <= value_threshold: + idx_ascend = idx + break + idx += 1 + + # If there is a single peak at fermi energy, no + # search needs to be performed. + if idx_ascend != fermi_idx and idx_descend != fermi_idx: + self.m_cache['highest_occupied_energy'] = fermi_energy_closest + self.m_cache['lowest_unoccupied_energy'] = fermi_energy_closest + single_peak_fermi = True + + if not single_peak_fermi: + # Look for highest occupied energy below the descend index + idx = idx_descend + while True: + try: + value = dos_values[idx] + except IndexError: + break + if value > value_threshold: + idx = idx if idx == idx_descend else idx + 1 + self.m_cache['highest_occupied_energy'] = energies[idx] + break + idx -= 1 + # Look for lowest unoccupied energy above idx_ascend + idx = idx_ascend + while True: + try: + value = dos_values[idx] + except IndexError: + break + if value > value_threshold: + idx = idx if idx == idx_ascend else idx - 1 + self.m_cache['highest_occupied_energy'] = energies[idx] + break + idx += 1 + + # Return the `highest_occupied_energy` as the `energies_origin`, or the `fermi_level` if it is not None + energies_origin = self.m_cache.get('highest_occupied_energy') + if energies_origin is None: + energies_origin = fermi_level + return energies_origin + + def resolve_normalization_factor(self, logger: BoundLogger) -> Optional[float]: + """ + Resolve the `normalization_factor` for the electronic DOS to get a cell-independent intensive DOS. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[float]): The normalization factor. + """ + # Get the `ModelSystem` as referenced in the `Outputs.model_system_ref` + model_system = get_sibling_section( + section=self, sibling_section_name='model_system_ref', logger=logger + ) + if model_system is None: + logger.warning( + 'Could not resolve the referenced `ModelSystem` in the `Outputs`.' + ) + return None + + # Get the originally parsed `AtomicCell`, which is the first element stored in `ModelSystem.cell` of name `'AtomicCell'` + atomic_cell = None + for cell in model_system.cell: + if cell.name == 'AtomicCell': # we get the originally parsed `AtomicCell` + atomic_cell = cell + break + if atomic_cell is None: + logger.warning( + 'Could not resolve the `AtomicCell` from the referenced `ModelSystem`.' + ) + return None + + # Get the `atoms_state` and their `atomic_number` from the `AtomicCell` + if atomic_cell.atoms_state is None or len(atomic_cell.atoms_state) == 0: + logger.warning('Could not resolve the `atoms_state` from the `AtomicCell`.') + return None + atomic_numbers = [atom.atomic_number for atom in atomic_cell.atoms_state] + + # Return `normalization_factor` depending if the calculation is spin polarized or not + if self.spin_channel is not None: + normalization_factor = 1 / (2 * sum(atomic_numbers)) + else: + normalization_factor = 1 / sum(atomic_numbers) + return normalization_factor + + def extract_band_gap(self) -> Optional[ElectronicBandGap]: + """ + Extract the electronic band gap from the `highest_occupied_energy` and `lowest_unoccupied_energy` stored + in `m_cache` from `resolve_energies_origin()`. If the difference of `highest_occupied_energy` and + `lowest_unoccupied_energy` is negative, the band gap `value` is set to 0.0. + + Returns: + (Optional[ElectronicBandGap]): The extracted electronic band gap section to be stored in `Outputs`. + """ + band_gap = None + homo = self.m_cache.get('highest_occupied_energy') + lumo = self.m_cache.get('lowest_unoccupied_energy') + if homo and lumo: + band_gap = ElectronicBandGap() + band_gap.is_derived = True + band_gap.physical_property_ref = self + + if (homo - lumo).magnitude < 0: + band_gap.value = 0.0 + else: + band_gap.value = homo - lumo + return band_gap + + def extract_projected_dos( + self, type: str, logger: BoundLogger + ) -> List[Optional[DOSProfile]]: + """ + Extract the projected DOS from the `projected_dos` section and the specified `type`. + + Args: + type (str): The type of the projected DOS to extract. It can be `'atom'` or `'orbital'`. + + Returns: + (DOSProfile): The extracted projected DOS. + """ + extracted_pdos = [] + for pdos in self.projected_dos: + # We make sure each PDOS is normalized + pdos.normalize(None, logger) + + # Initial check for `name` and `entity_ref` + if ( + pdos.name is None + or pdos.entity_ref is None + or len(pdos.entity_ref) == 0 + ): + logger.warning( + '`name` or `entity_ref` are not set for `projected_dos` and they are required for normalization to work.' + ) + return None + + if type in pdos.name: + extracted_pdos.append(pdos) + return extracted_pdos + + def generate_from_projected_dos( + self, logger: BoundLogger + ) -> Optional[pint.Quantity]: + """ + Generate the total `value` of the electronic DOS from the `projected_dos` contributions. If the `projected_dos` + is not present, it returns `None`. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[pint.Quantity]): The total `value` of the electronic DOS. + """ + if self.projected_dos is None or len(self.projected_dos) == 0: + return None + + # Extract `Energy` variables + energies = get_variables(self.variables, Energy) + if len(energies) != 1: + logger.warning( + 'The `ElectronicDensityOfStates` does not contain an `Energy` variable to extract the DOS.' + ) + return None + + # We distinguish between orbital and atom `projected_dos` + orbital_projected = self.extract_projected_dos('orbital', logger) + atom_projected = self.extract_projected_dos('atom', logger) + + # Extract `atom_projected` from `orbital_projected` by summing up the `orbital_projected` contributions for each atom + if len(atom_projected) == 0: + atom_data: Dict[AtomsState, List[DOSProfile]] = {} + for orb_pdos in orbital_projected: + # `entity_ref` is the `OrbitalsState` section, whose parent is `AtomsState` + entity_ref = orb_pdos.entity_ref.m_parent + if entity_ref in atom_data: + atom_data[entity_ref].append(orb_pdos) + else: + atom_data[entity_ref] = [orb_pdos] + for ref, data in atom_data.items(): + atom_dos = DOSProfile( + name=f'atom {ref.chemical_symbol}', + entity_ref=ref, + variables=energies, + ) + orbital_values = [ + dos.value.magnitude for dos in data + ] # to avoid warnings from pint + orbital_unit = data[0].value.u + atom_dos.value = np.sum(orbital_values, axis=0) * orbital_unit + atom_projected.append(atom_dos) + # We concatenate the `atom_projected` to the `projected_dos` + self.projected_dos = orbital_projected + atom_projected + + # Extract `value` from `atom_projected` by summing up the `atom_projected` contributions + value = self.value + if value is None: + atom_values = [ + dos.value.magnitude for dos in atom_projected + ] # to avoid warnings from pint + atom_unit = atom_projected[0].value.u + value = np.sum(atom_values, axis=0) * atom_unit + return value + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Initial check to see if `variables` contains the required `Energy` variable + energies = get_variables(self.variables, Energy) + if len(energies) != 0: + return + energies = energies[0].points + + # Resolve `fermi_level` from a sibling section with respect to `ElectronicDensityOfStates` + fermi_level = get_sibling_section( + section=self, sibling_section_name='fermi_level', logger=logger + ) # * we consider `index_sibling` to be 0 + if fermi_level is not None: + fermi_level = fermi_level.value + # and the `energies_origin` from the sibling `ElectronicEigenvalues` section + self.energies_origin = self.resolve_energies_origin( + energies, fermi_level, logger + ) + if self.energies_origin is None: + logger.info('Could not resolve the `energies_origin` for the DOS') + + # Resolve `normalization_factor` + if self.normalization_factor is None: + self.normalization_factor = self.resolve_normalization_factor(logger) + + # `ElectronicBandGap` extraction + band_gap = self.extract_band_gap() + if band_gap is not None: + self.m_parent.electronic_band_gap.append(band_gap) + + # Total `value` extraction from `projected_dos` + value_from_pdos = self.generate_from_projected_dos(logger) + if self.value is None and value_from_pdos is not None: + logger.info( + 'The `ElectronicDensityOfStates.value` is not stored. We will attempt to obtain it by summing up projected DOS contributions, if these are present.' + ) + self.value = value_from_pdos + + +class AbsorptionSpectrum(SpectralProfile): + """ """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + axis = Quantity( + type=MEnum('xx', 'yy', 'zz'), + description=""" + Axis of the absorption spectrum. This is related with the polarization direction, and can be seen as the + principal term in the tensor `Permittivity.value` (see permittivity.py module). + """, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # Set the name of the section + self.name = self.m_def.name + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class XASSpectrum(AbsorptionSpectrum): + """ + X-ray Absorption Spectrum (XAS). + """ + + # ! implement `iri` and `rank` as part of `m_def = Section()` + + xanes_spectrum = SubSection( + sub_section=AbsorptionSpectrum.m_def, + description=""" + X-ray Absorption Near Edge Structure (XANES) spectrum. + """, + repeats=False, + ) + + exafs_spectrum = SubSection( + sub_section=AbsorptionSpectrum.m_def, + description=""" + Extended X-ray Absorption Fine Structure (EXAFS) spectrum. + """, + repeats=False, + ) + + def __init__( + self, m_def: Section = None, m_context: Context = None, **kwargs + ) -> None: + super().__init__(m_def, m_context, **kwargs) + # Set the name of the section + self.name = self.m_def.name + + def generate_from_contributions(self, logger: BoundLogger) -> None: + """ + Generate the `value` of the XAS spectrum by concatenating the XANES and EXAFS contributions. It also concatenates + the `Energy` grid points of the XANES and EXAFS parts. + + Args: + logger (BoundLogger): The logger to log messages. + """ + # TODO check if this method is general enough + if self.xanes_spectrum is not None and self.exafs_spectrum is not None: + # Concatenate XANE and EXAFS `Energy` grid points + xanes_variables = get_variables(self.xanes_spectrum.variables, Energy) + exafs_variables = get_variables(self.exafs_spectrum.variables, Energy) + if len(xanes_variables) == 0 or len(exafs_variables) == 0: + logger.warning( + 'Could not extract the `Energy` grid points from XANES or EXAFS.' + ) + return + xanes_energies = xanes_variables[0].points + exafs_energies = exafs_variables[0].points + if xanes_energies.max() > exafs_energies.min(): + logger.warning( + 'The XANES `Energy` grid points are not below the EXAFS `Energy` grid points.' + ) + return + self.variables = [ + Energy(points=np.concatenate([xanes_energies, exafs_energies])) + ] + # Concatenate XANES and EXAFS `value` if they have the same shape ['n_energies'] + try: + self.value = np.concatenate( + [self.xanes_spectrum.value, self.exafs_spectrum.value] + ) + except ValueError: + logger.warning( + 'The XANES and EXAFS `value` have different shapes. Could not concatenate the values.' + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + if self.value is None: + logger.info( + 'The `XASSpectrum.value` is not stored. We will attempt to obtain it by combining the XANES and EXAFS parts if these are present.' + ) + self.generate_from_contributions(logger) + diff --git a/properties/thermodynamics.py b/properties/thermodynamics.py new file mode 100644 index 00000000..fdd05883 --- /dev/null +++ b/properties/thermodynamics.py @@ -0,0 +1,239 @@ +# +# 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. +# +# +# 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, SubSection, MEnum +from .physical_property import PhysicalProperty + +class Enthalpy(PhysicalProperty): + """ + Section containing the enthalpy (i.e. energy_total + pressure * volume.) of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the enthalpy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Entropy(PhysicalProperty): + """ + Section containing the entropy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule / kelvin', + description=""" + Value of the entropy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ChemicalPotential(PhysicalProperty): + """ + Section containing the chemical potential of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the chemical potential. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Pressure(PhysicalProperty): + """ + Section containing the pressure of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='pascal', + description=""" + Value of the pressure. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Virial(PhysicalProperty): + """ + Section containing the virial (cross product between positions and forces) of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the pressure. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Temperature(PhysicalProperty): + """ + Section containing the temperature of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='kelvin', + description=""" + Value of the pressure. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Volume(PhysicalProperty): + """ + Section containing the volume of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='m ** 3', + description=""" + Value of the volume. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Density(PhysicalProperty): + """ + Section containing the density of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='kg / m ** 3', + description=""" + Value of the density. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +# ? Does this go here or in energies? +# ? Naming specific to Potential Energy? +class Hessian(PhysicalProperty): + """ + Section containing the Hessian matrix, i.e., 2nd derivatives with respect to geometric (typically particle) displacements, + of the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule / m ** 2', + description=""" + Value of the Hessian. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class HeatCapacityCV(PhysicalProperty): + """ + Section containing the heat capacity at constant volume for a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule / kelvin', + description=""" + Value of the heat capacity. + """, + ) + +class HeatCapacityCP(PhysicalProperty): + """ + Section containing the heat capacity at constant volume for a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule / kelvin', + description=""" + Value of the heat capacity. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + + # ? Is this ever used? + # internal_energy = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Value of the internal energy. + # """, + # ) + + # vibrational_free_energy_at_constant_volume = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Value of the vibrational free energy per cell unit at constant volume. + # """, + # ) + + + diff --git a/src/nomad_simulations/outputs.py b/src/nomad_simulations/outputs.py index 388e6504..dae9e87e 100644 --- a/src/nomad_simulations/outputs.py +++ b/src/nomad_simulations/outputs.py @@ -24,6 +24,7 @@ from nomad.datamodel.metainfo.annotations import ELNAnnotation from nomad_simulations.model_system import ModelSystem +from nomad_simulations.model_method import ModelMethod from nomad_simulations.physical_property import PhysicalProperty from nomad_simulations.numerical_settings import SelfConsistency from nomad_simulations.properties import ( @@ -39,6 +40,7 @@ AbsorptionSpectrum, XASSpectrum, Permittivity, + TotalEnergy ) @@ -64,6 +66,15 @@ class Outputs(ArchiveSection): a_eln=ELNAnnotation(component='ReferenceEditQuantity'), ) + model_method_ref = Quantity( + type=ModelMethod, + description=""" + Reference to the `ModelMethod` section to which the output property references to + and on which the simulation is performed. + """, + a_eln=ELNAnnotation(component='ReferenceEditQuantity'), + ) + # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # List of properties # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # @@ -100,6 +111,8 @@ class Outputs(ArchiveSection): xas_spectra = SubSection(sub_section=XASSpectrum.m_def, repeats=True) + total_energy = SubSection(SubSection=TotalEnergy.m_def, repeats=True) + # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # @@ -281,3 +294,43 @@ def normalize(self, archive, logger) -> None: physical_property=phys_property, logger=logger, ) + +class WorkflowOutputs(Outputs): + """ + This section contains output properties that depend on a single system, but were + calculated as part of a workflow. Examples include geometry optimization and molecular dynamics. + """ + + step = Quantity( + type=np.int32, + description=""" + The step number with respect to the workflow. + """, + ) + + workflow_ref = Quantity( + type=SimulationWorkflow, + description=""" + Reference to the `SelfConsistency` section that defines the numerical settings to converge the + output property. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class TrajectoryOutputs(WorkflowOutputs): + """ + This section contains output properties that depend on a single system, but were + calculated as part of a workflow. Examples include geometry optimization and molecular dynamics. + """ + + time = Quantity( + type=np.float64, + description=""" + The elapsed simulated physical time since the start of the trajectory. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) diff --git a/src/nomad_simulations/physical_property.py b/src/nomad_simulations/physical_property.py index 39c1d926..e879a593 100644 --- a/src/nomad_simulations/physical_property.py +++ b/src/nomad_simulations/physical_property.py @@ -320,4 +320,4 @@ def normalize(self, archive, logger) -> None: super().normalize(archive, logger) # Resolve if the physical property `is_derived` or not from another physical property. - self.is_derived = self._is_derived() + self.is_derived = self._is_derived() \ No newline at end of file diff --git a/src/nomad_simulations/properties/__init__.py b/src/nomad_simulations/properties/__init__.py index 3c6194ed..0baae70b 100644 --- a/src/nomad_simulations/properties/__init__.py +++ b/src/nomad_simulations/properties/__init__.py @@ -16,7 +16,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -from .energies import FermiLevel, ChemicalPotential +from .energies import FermiLevel, ChemicalPotential, TotalEnergy from .band_gap import ElectronicBandGap from .spectral_profile import ( SpectralProfile, diff --git a/src/nomad_simulations/properties/energies.py b/src/nomad_simulations/properties/energies.py index 3283a1f8..b728c97c 100644 --- a/src/nomad_simulations/properties/energies.py +++ b/src/nomad_simulations/properties/energies.py @@ -18,11 +18,10 @@ import numpy as np -from nomad.metainfo import Quantity, Section, Context +from nomad.metainfo import Quantity, Section, Context, SubSection, MEnum from nomad_simulations.physical_property import PhysicalProperty - class FermiLevel(PhysicalProperty): """ Energy required to add or extract a charge from a material at zero temperature. It can be also defined as the chemical potential at zero temperature. @@ -77,3 +76,675 @@ def __init__( def normalize(self, archive, logger) -> None: super().normalize(archive, logger) + + +#################################################### +# List of classical energy contribuions +#################################################### + +class PotentialEnergy(PhysicalProperty): + """ + Section containing the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class KineticEnergy(PhysicalProperty): + """ + Section containing the kinetic energy of a (sub)system. + """ + + type = Quantity( + type=MEnum('classical', 'quantum'), + description=""" + Refers to the method used for calculating the kinetic energy. + + Allowed values are: + + | Thermostat Name | Description | + + | ---------------------- | ----------------------------------------- | + + | `"classical"` | The kinetic energy is calculated directly from + the velocities of particles as KE = /sum_i 1/2 m_i v_i^2, where the sum runs over the number of particles in the system, + m_i is the mass of particle i, and v_i is the velocity of particle i. | + + | `"quantum"` | ... | + """, + ) + # ? Just an idea, not sure if this is necessary since we will be referencing the method, I guess maybe it's overkill in general + # ? If we do this, maybe we should have a general Energy(PP) class which is inherited from in case there are multiple interpretations/calculation methods + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the kinetic energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class VDWEnergy(PhysicalProperty): + """ + Section containing the van der Waals contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the van der Waals contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectrostaticEnergy(PhysicalProperty): + """ + Section containing all electrostatic contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of all electrostatic contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectrostaticShortRangeEnergy(PhysicalProperty): + """ + Section containing short-range electrostatic contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the short-range electrostatic contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectrostaticLongRangeEnergy(PhysicalProperty): + """ + Section containing long-range electrostatic contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the long-range electrostatic contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class BondedEnergy(PhysicalProperty): + """ + Section containing all bonded (i.e., intramolecular) contributions to the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of all bonded (i.e., intramolecular) contributions to the potential energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class BondEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from bond interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from bond interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class AngleEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from angle interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from angle interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class DihedralEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from dihedral interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from dihedral interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ImproperDihedralEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from improper dihedral interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from improper dihedral interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ExternalEnergy(PhysicalProperty): + """ + Section containing contributions to the potential energy from external interactions of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of contributions to the potential energy from external interactions. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ClassicalEnergyContributions(PhysicalProperty): + """ + Section containing contributions to the potential energy from a classical force field. + """ + + potential = SubSection(sub_section=PotentialEnergy.m_def, repeats=False) + + kinetic = SubSection(sub_section=KineticEnergy.m_def, repeats=False) + + vdw = SubSection(sub_section=VDWEnergy.m_def, repeats=False) + + electrostatic = SubSection(sub_section=ElectrostaticEnergy.m_def, repeats=False) + + electrostatic_short_range = SubSection(sub_section=ElectrostaticShortRangeEnergy.m_def, repeats=False) + + electrostatic_long_range = SubSection(sub_section=ElectrostaticLongRangeEnergy.m_def, repeats=False) + + bonded = SubSection(sub_section=BondedEnergy.m_def, repeats=False) + + bond = SubSection(sub_section=BondEnergy.m_def, repeats=False) + + angle = SubSection(sub_section=AngleEnergy.m_def, repeats=False) + + dihedral = SubSection(sub_section=DihedralEnergy.m_def, repeats=False) + + improper_dihedral = SubSection(sub_section=ImproperDihedralEnergy.m_def, repeats=False) + + external = SubSection(sub_section=ExternalEnergy.m_def, repeats=False) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + # Set the name of the section + self.name = self.m_def.name + +###################################### +# List of quantum energy contributions +###################################### + +class ElectronicEnergy(PhysicalProperty): + """ + Section containing the electronic energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the electronic energy. + """, + ) + + contributions = SubSection(sub_section=EnergyContributions.m_def, repeats=True) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + allowed_contributions = ['PotKin'] + +class ElectronicKineticEnergy(PhysicalProperty): + """ + Section containing the electronic kinetic energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the electronic kinetic energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class XCEnergy(PhysicalProperty): + """ + Section containing the exchange-correlation (XC) energy of a (sub)system, + calculated using the functional stored in XC_functional. + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the exchange-correlation energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class XCPotentialEnergy(PhysicalProperty): + """ + Section containing the potential energy contribution to the exchange-correlation (XC) energy, + i.e., the integral of the first order derivative of the functional + stored in XC_functional (integral of v_xc*electron_density), i.e., the component + of XC that is in the sum of the eigenvalues. Value associated with the + configuration, should be the most converged value. | + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the potential energy contribution of the exchange-correlation energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +# ? XCCorrelationEnergy? +class CorrelationEnergy(PhysicalProperty): + """ + Section containing the correlation energy of a (sub)system, + calculated using the method described in XC_functional. + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the correlation energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +# ? XCExchangeEnergy? +class ExchangeEnergy(PhysicalProperty): + """ + Section containing the exchange energy of a (sub)system, + calculated using the method described in XC_functional. + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the exchange energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class T0Energy(PhysicalProperty): + """ + Section containing the total energy of a (sub)system extrapolated to $T=0$. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the total energy extrapolated to + $T=0$, based on a free-electron gas argument. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ZeroPointEnergy(PhysicalProperty): + """ + Section containing the zero-point vibrational energy of a (sub)system, + calculated using the method described in zero_point_method. + """ + # ! Someone check this description! + # ? Do we really want to specify the method here? This can't be user-defined? + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the zero-point energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ElectrostaticEnergy(PhysicalProperty): + """ + Section containing the electrostatic energy (nuclei + electrons) of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the electrostatic energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class NuclearRepulsionEnergy(PhysicalProperty): + """ + Section containing the nuclear-nuclear repulsion energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the nuclear-nuclear repulsion energy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class QuantumEnergyContributions(EnergyContributions): + """ + Section containing contributions to the potential energy from a DFT calculation. + """ +... + + electronic = SubSection(sub_section=ElectronicEnergy.m_def, repeats=False) + + electronic_kinetic = SubSection(sub_section=ElectronicKineticEnergy.m_def, repeats=False) + + xc = SubSection(sub_section=XCEnergy.m_def, repeats=False) + + xc_potential = SubSection(sub_section=XCPotentialEnergy.m_def, repeats=False) + + correlation = SubSection(sub_section=CorrelationEnergy.m_def, repeats=False) + + exchange = SubSection(sub_section=ExchangeEnergy.m_def, repeats=False) + + t0 = SubSection(sub_section=T0Energy.m_def, repeats=False) + + zero_point = SubSection(sub_section=ZeroPointEnergy.m_def, repeats=False) + + electrostatic = SubSection(sub_section=ElectrostaticEnergy.m_def, repeats=False) + + nuclear_repulsion = SubSection(sub_section=NuclearRepulsionEnergy.m_def, repeats=False) + +########################## +# Other / General energies +########################## + +class TotalEnergy(PhysicalProperty): + """ + Section containing the total energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + The value of the total energy. + """, + ) + # ? Do we need these descriptions under value? It ends up simply duplicating the section info to some extent. + + classical_contributions = SubSection(sub_section=ClassicalEnergyContributions.m_def, repeats=False) + + quantum_contributions = SubSection(sub_section=QuantumEnergyContributions.m_def, repeats=False) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # madelung = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the Madelung energy. + # """, + # ) + + # # TODO I suggest ewald is moved to "long range" under electrostatic->energyentry, unless there is some other usage I am misunderstanding + # ewald = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the Ewald energy. + # """, + # ) + + # free = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the free energy (nuclei + electrons) + # (whose minimum gives the smeared occupation density calculated with + # smearing_kind). + # """, + # ) + + # sum_eigenvalues = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the sum of the eigenvalues of the + # Hamiltonian matrix. + # """, + # ) + + # van_der_waals = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the Van der Waals energy. A multiple + # occurence is expected when more than one van der Waals methods are defined. The + # van der Waals kind should be specified in Energy.kind + # """, + # ) + + # hartree_fock_x_scaled = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Scaled exact-exchange energy that depends on the mixing parameter of the + # functional. For example in hybrid functionals, the exchange energy is given as a + # linear combination of exact-energy and exchange energy of an approximate DFT + # functional; the exact exchange energy multiplied by the mixing coefficient of the + # hybrid functional would be stored in this metadata. Defined consistently with + # XC_method. + # """, + # ) + + # ? This is technically NOT redundant with total energy, but I fear that people will use them interchangeably, so we need to be clear about the definitions in any case + # internal = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Value of the internal energy. + # """, + # ) + + # double_counting = SubSection( + # sub_section=EnergyEntry.m_def, + # categories=[FastAccess], + # description=""" + # Double counting correction when performing Hubbard model calculations. + # """, + # ) + + # # TODO remove this should be be entropy.correction + # correction_entropy = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Entropy correction to the potential energy to compensate for the change in + # occupation so that forces at finite T do not need to keep the change of occupation + # in account. Defined consistently with XC_method. + # """, + # ) + + # # TODO remove this should be in electrostatic.correction + # correction_hartree = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Correction to the density-density electrostatic energy in the sum of eigenvalues + # (that uses the mixed density on one side), and the fully consistent density- + # density electrostatic energy. Defined consistently with XC_method. + # """, + # ) + + # # TODO remove this should be in xc.correction + # correction_xc = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Correction to energy_XC. + # """, + # ) + + # # ? Is it ok to store this in outputs and not in workflow? I guess we can ensure in normalization that this is a WorkflowOutput, etc...? + # change = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Stores the change of total energy with respect to the previous step. + # """, + # categories=[ErrorEstimateContribution, EnergyValue], + # ) + + # fermi = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Fermi energy (separates occupied from unoccupied single-particle states) + # """, + # categories=[EnergyTypeReference, EnergyValue], + # ) + + # highest_occupied = Quantity( + # type=np.dtype(np.float64), + # unit='joule', + # shape=[], + # description=""" + # The highest occupied energy. + # """, + # ) + + # lowest_unoccupied = Quantity( + # type=np.dtype(np.float64), + # unit='joule', + # shape=[], + # description=""" + # The lowest unoccupied energy. + # """, + # ) + + + # # TODO this should be removed and replaced by correction in EnergyEntry + # current = SubSection( + # sub_section=EnergyEntry.m_def, + # description=""" + # Contains the value and information regarding the energy calculated with + # calculation_method_current. energy_current is equal to energy_total for + # non-perturbative methods. For perturbative methods, energy_current is equal to the + # correction: energy_total minus energy_total of the calculation_to_calculation_ref + # with calculation_to_calculation_kind = starting_point + # """, + # ) + + + + +# # ? Do we want to allow this? +# class EnergyCustom(PhysicalProperty): +# """ +# Section containing the total energy of a (sub)system. +# """ + +# name = Quantity( +# type=str, +# description=""" +# The name of this custom energy. +# """, +# ) + +# value = Quantity( +# type=np.float64, +# unit='joule', +# description=""" +# The value of this custom energy type. +# """, +# ) + +# def normalize(self, archive, logger) -> None: +# super().normalize(archive, logger) diff --git a/src/nomad_simulations/properties/forces.py b/src/nomad_simulations/properties/forces.py new file mode 100644 index 00000000..8600f190 --- /dev/null +++ b/src/nomad_simulations/properties/forces.py @@ -0,0 +1,181 @@ +# +# 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. +# +# +# 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 re +import numpy as np + +from nomad.units import ureg +from nomad.datamodel.data import ArchiveSection +from nomad.metainfo import Quantity, SubSection, SectionProxy, MEnum +from nomad.metainfo.metainfo import DirectQuantity, Dimension +from nomad.datamodel.metainfo.annotations import ELNAnnotation + +from .outputs import Outputs, SCFOutputs, WorkflowOutputs, TrajectoryOutputs +from .atoms_state import AtomsState +from .physical_property import PhysicalProperty + + +from nomad.datamodel.metainfo.common import ( + FastAccess, + PropertySection, + ProvenanceTracker, +) +from nomad.metainfo import ( + Category, + HDF5Reference, + MCategory, + MEnum, + MSection, + Package, + Quantity, + Reference, + Section, + SectionProxy, + SubSection, +) + +from .model_system import ModelSystem + + + +class ForceTotal(PhysicalProperty): + """ + Section containing the total force of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='newton', + description=""" + The value of the total force. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + + +class ForcesEntry(Atomic): + """ + Section describing a contribution to or type of atomic forces. + """ + + m_def = Section(validate=False) + + value = Quantity( + type=np.dtype(np.float64), + shape=['n_atoms', 3], + unit='newton', + description=""" + Value of the forces acting on the atoms. This is calculated as minus gradient of + the corresponding energy type or contribution **including** constraints, if + present. The derivatives with respect to displacements of nuclei are evaluated in + Cartesian coordinates. In addition, these are obtained by filtering out the + unitary transformations (center-of-mass translations and rigid rotations for + non-periodic systems, see value_raw for the unfiltered counterpart). + """, + ) + + value_raw = Quantity( + type=np.dtype(np.float64), + shape=['n_atoms', 3], + unit='newton', + description=""" + Value of the forces acting on the atoms **not including** such as fixed atoms, + distances, angles, dihedrals, etc.""", + ) + # ? This is VERY imprecise, is this used regularly? + +class Forces(MSection): + """ + Section containing all forces types and contributions. + """ + + m_def = Section(validate=False) + + total = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains the value and information regarding the total forces on the atoms + calculated as minus gradient of energy_total. + """, + ) + # ! We need to avoid giving the precise method of calculation without also providing context, this is not necessarily true in general! + + free = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains the value and information regarding the forces on the atoms + corresponding to the minus gradient of energy_free. The (electronic) energy_free + contains the information on the change in (fractional) occupation of the + electronic eigenstates, which are accounted for in the derivatives, yielding a + truly energy-conserved quantity. + """, + ) + + t0 = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains the value and information regarding the forces on the atoms + corresponding to the minus gradient of energy_T0. + """, + ) + + contributions = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains other forces contributions to the total atomic forces not already + defined. + """, + repeats=True, + ) + + types = SubSection( + sub_section=ForcesEntry.m_def, + description=""" + Contains other types of forces not already defined. + """, + repeats=True, + ) + + + + + + forces = SubSection(sub_section=Forces.m_def) + + diff --git a/src/nomad_simulations/properties/mechanics.py b/src/nomad_simulations/properties/mechanics.py new file mode 100644 index 00000000..bbe418ba --- /dev/null +++ b/src/nomad_simulations/properties/mechanics.py @@ -0,0 +1,129 @@ +# +# 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. +# +# +# 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 re +import numpy as np + +from nomad.units import ureg +from nomad.datamodel.data import ArchiveSection +from nomad.metainfo import Quantity, SubSection, SectionProxy, MEnum +from nomad.metainfo.metainfo import DirectQuantity, Dimension +from nomad.datamodel.metainfo.annotations import ELNAnnotation + +from .outputs import Outputs, SCFOutputs, WorkflowOutputs, TrajectoryOutputs +from .atoms_state import AtomsState +from .physical_property import PhysicalProperty + + +from nomad.datamodel.metainfo.common import ( + FastAccess, + PropertySection, + ProvenanceTracker, +) +from nomad.metainfo import ( + Category, + HDF5Reference, + MCategory, + MEnum, + MSection, + Package, + Quantity, + Reference, + Section, + SectionProxy, + SubSection, +) + +from .model_system import ModelSystem + + +class StressEntry(Atomic): + """ + Section describing a contribution to or a type of stress. + """ + + m_def = Section(validate=False) + + value = Quantity( + type=np.dtype(np.float64), + shape=[3, 3], + unit='joule/meter**3', + description=""" + Value of the stress on the simulation cell. It is given as the functional + derivative of the corresponding energy with respect to the deformation tensor. + """, + ) + + values_per_atom = Quantity( + type=np.dtype(np.float64), + shape=['number_of_atoms', 3, 3], + unit='joule/meter**3', + description=""" + Value of the atom-resolved stresses. + """, + ) + + +class Stress(MSection): + """ + Section containing all stress types and contributions. + """ + + m_def = Section(validate=False) + + total = SubSection( + sub_section=StressEntry.m_def, + description=""" + Contains the value and information regarding the stress on the simulation cell + and the atomic stresses corresponding to energy_total. + """, + ) + + contributions = SubSection( + sub_section=StressEntry.m_def, + description=""" + Contains contributions for the total stress. + """, + repeats=True, + ) + + types = SubSection( + sub_section=StressEntry.m_def, + description=""" + Contains other types of stress. + """, + repeats=True, + ) \ No newline at end of file diff --git a/src/nomad_simulations/properties/molecular.py b/src/nomad_simulations/properties/molecular.py new file mode 100644 index 00000000..931e0930 --- /dev/null +++ b/src/nomad_simulations/properties/molecular.py @@ -0,0 +1,109 @@ +# +# 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. +# +# +# 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 re +import numpy as np + +from nomad.units import ureg +from nomad.datamodel.data import ArchiveSection +from nomad.metainfo import Quantity, SubSection, SectionProxy, MEnum +from nomad.metainfo.metainfo import DirectQuantity, Dimension +from nomad.datamodel.metainfo.annotations import ELNAnnotation + +from .outputs import Outputs, SCFOutputs, WorkflowOutputs, TrajectoryOutputs +from .atoms_state import AtomsState +from .physical_property import PhysicalProperty + + +class RadiusOfGyration(PhysicalProperty): + + value = Quantity( + type=np.float64, + unit='m', + ) + +class Hessian(TrajectoryOutputs): + """ + Section containing the Hessian matrix, i.e., 2nd derivatives with respect to atomic displacements, + of the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.dtype(np.float64), + shape=['n_atoms', 'n_atoms', 3, 3], + unit='joule / m ** 2', + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + self.value_unit = 'joule / m ** 2' + self.name = 'hessian' + self.variables = [['atom_index'], ['atom_index'], ['x', 'y', 'z'], ['x', 'y', 'z']] + self.bins = [np.range(self.n_atoms), np.range(self.n_atoms)] + + +class RadialDistributionFunction(TrajectoryOutputs): + """ + Section containing information about the calculation of + radial distribution functions (rdfs). + """ + + value = Quantity( + type=np.dtype(np.float64), + shape=['*'], + ) + + variables = Quantity( + type=str, + shape=[1], + description=""" + Name/description of the variables along which the property is defined. + """, + ) + + bins = Quantity( + type=np.float64, + shape=['*'], + unit='m', + description=""" + Distances along which the rdf was calculated. + """, + ) + + def normalize(self, archive, logger): + super().normalize(archive, logger) + self.variables = ['distance'] + assert len(self.bins) == len(self.value) diff --git a/src/nomad_simulations/properties/thermodynamics.py b/src/nomad_simulations/properties/thermodynamics.py new file mode 100644 index 00000000..ecd1ec62 --- /dev/null +++ b/src/nomad_simulations/properties/thermodynamics.py @@ -0,0 +1,237 @@ +# +# 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. +# +# +# 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, SubSection, MEnum +from .physical_property import PhysicalProperty + +class Enthalpy(PhysicalProperty): + """ + Section containing the enthalpy (i.e. energy_total + pressure * volume.) of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the enthalpy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Entropy(PhysicalProperty): + """ + Section containing the entropy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule / kelvin', + description=""" + Value of the entropy. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ChemicalPotential(PhysicalProperty): + """ + Section containing the chemical potential of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the chemical potential. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Pressure(PhysicalProperty): + """ + Section containing the pressure of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='pascal', + description=""" + Value of the pressure. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Virial(PhysicalProperty): + """ + Section containing the virial (cross product between positions and forces) of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule', + description=""" + Value of the pressure. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Temperature(PhysicalProperty): + """ + Section containing the temperature of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='kelvin', + description=""" + Value of the pressure. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Volume(PhysicalProperty): + """ + Section containing the volume of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='m ** 3', + description=""" + Value of the volume. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Density(PhysicalProperty): + """ + Section containing the density of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='kg / m ** 3', + description=""" + Value of the density. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +# ? Does this go here or in energies? +# ? Naming specific to Potential Energy? +class Hessian(PhysicalProperty): + """ + Section containing the Hessian matrix, i.e., 2nd derivatives with respect to geometric (typically particle) displacements, + of the potential energy of a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule / m ** 2', + description=""" + Value of the Hessian. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class HeatCapacityCV(PhysicalProperty): + """ + Section containing the heat capacity at constant volume for a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule / kelvin', + description=""" + Value of the heat capacity. + """, + ) + +class HeatCapacityCP(PhysicalProperty): + """ + Section containing the heat capacity at constant volume for a (sub)system. + """ + + value = Quantity( + type=np.float64, + unit='joule / kelvin', + description=""" + Value of the heat capacity. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + + # ? Is this ever used? + # internal_energy = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Value of the internal energy. + # """, + # ) + + # vibrational_free_energy_at_constant_volume = Quantity( + # type=np.dtype(np.float64), + # shape=[], + # unit='joule', + # description=""" + # Value of the vibrational free energy per cell unit at constant volume. + # """, + # ) + diff --git a/src/nomad_simulations/test_direct.py b/src/nomad_simulations/test_direct.py new file mode 100644 index 00000000..88eec4cf --- /dev/null +++ b/src/nomad_simulations/test_direct.py @@ -0,0 +1,48 @@ +# +# 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.units import ureg + +from nomad_simulations.outputs import ElectronicBandGap + +from nomad_simulations.variables import Variables, Temperature # ? delete these imports + +# Playing with `PhysicalProperty` +band_gap = ElectronicBandGap(source='simulation', type='direct', label='DFT') +n_grid_points = 3 +temperature = Temperature(n_grid_points=n_grid_points, grid_points=np.linspace(0, 100, n_grid_points)) +band_gap.variables.append(temperature) +n_grid_points = 2 +custom_bins = Variables(n_grid_points=n_grid_points, grid_points=np.linspace(0, 100, n_grid_points)) +band_gap.variables.append(custom_bins) +# band_gap.value_unit = 'joule' +# band_gap.shape = [3] +# band_gap.value = [ +# [[1, 2, 3], [1, 2, 3]], +# [[1, 2, 3], [1, 2, 3]], +# [[1, 2, 3], [1, 2, 3]], +# ] * ureg.eV +band_gap.value = [ + [1, 1], + [1, 1], + [1, 1], +] * ureg.eV +# band_gap.value = [1, 2, 3] * ureg.eV +print(band_gap)