diff --git a/simulationdataschema/atoms_state.py b/simulationdataschema/atoms_state.py new file mode 100644 index 00000000..97ce0bcd --- /dev/null +++ b/simulationdataschema/atoms_state.py @@ -0,0 +1,598 @@ +# +# 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 +import ase +from typing import Optional +from structlog.stdlib import BoundLogger + +from nomad.units import ureg + +from nomad.metainfo import Quantity, SubSection, MEnum +from nomad.datamodel.data import ArchiveSection +from nomad.datamodel.metainfo.annotations import ELNAnnotation + +from .utils import RussellSaundersState + + +class OrbitalsState(ArchiveSection): + """ + A base section used to define the orbital state of an atom. + """ + + # TODO: add the relativistic kappa_quantum_number + n_quantum_number = Quantity( + type=np.int32, + description=""" + Principal quantum number of the orbital state. + """, + ) + + l_quantum_number = Quantity( + type=np.int32, + description=""" + Azimuthal quantum number of the orbital state. This quantity is equivalent to `l_quantum_symbol`. + """, + ) + + l_quantum_symbol = Quantity( + type=str, + description=""" + Azimuthal quantum symbol of the orbital state, "s", "p", "d", "f", etc. This + quantity is equivalent to `l_quantum_number`. + """, + ) + + ml_quantum_number = Quantity( + type=np.int32, + description=""" + Azimuthal projection number of the `l` vector. This quantity is equivalent to `ml_quantum_symbol`. + """, + ) + + ml_quantum_symbol = Quantity( + type=str, + description=""" + Azimuthal projection symbol of the `l` vector, "x", "y", "z", etc. This quantity is equivalent + to `ml_quantum_number`. + """, + ) + + j_quantum_number = Quantity( + type=np.float64, + shape=["1..2"], + description=""" + Total angular momentum quantum number $j = |l-s| ... l+s$. Necessary with strong + L-S coupling or non-collinear spin systems. + """, + ) + + mj_quantum_number = Quantity( + type=np.float64, + shape=["*"], + description=""" + Azimuthal projection of the `j` vector. Necessary with strong L-S coupling or + non-collinear spin systems. + """, + ) + + ms_quantum_number = Quantity( + type=np.int32, + description=""" + Spin quantum number. Set to -1 for spin down and +1 for spin up. In non-collinear spin + systems, the projection axis $z$ should also be defined. + """, + ) + + ms_quantum_symbol = Quantity( + type=str, + description=""" + Spin quantum symbol. Set to 'down' for spin down and 'up' for spin up. In non-collinear + spin systems, the projection axis $z$ should also be defined. + """, + ) + + degeneracy = Quantity( + type=np.int32, + description=""" + The degeneracy of the orbital state. The degeneracy is the number of states with the same + energy. It is equal to 2 * l + 1 for non-relativistic systems and 2 * j + 1 for + relativistic systems, if ms_quantum_number is defined (otherwise a factor of 2 is included). + """, + ) + + occupation = Quantity( + type=np.float64, + description=""" + The number of electrons in the orbital state. The state is fully occupied if + occupation = degeneracy. + """, + ) + + def __init__(self): + self._orbitals = { + -1: dict(zip(range(4), ("s", "p", "d", "f"))), + 0: {0: ""}, + 1: dict(zip(range(-1, 2), ("x", "z", "y"))), + 2: dict(zip(range(-2, 3), ("xy", "xz", "z^2", "yz", "x^2-y^2"))), + 3: dict( + zip( + range(-3, 4), + ( + "x(x^2-3y^2)", + "xyz", + "xz^2", + "z^3", + "yz^2", + "z(x^2-y^2)", + "y(3x^2-y^2)", + ), + ) + ), + } + self._orbitals_map = { + "l_symbols": self._orbitals[-1], + "ml_symbols": {i: self._orbitals[i] for i in range(4)}, + "ms_symbols": dict((zip((False, True), ("down", "up")))), + "l_numbers": {v: k for k, v in self._orbitals[-1].items()}, + "ml_numbers": { + k: {v: k for k, v in self._orbitals[k].items()} for k in range(4) + }, + "ms_numbers": dict((zip((False, True), (-0.5, 0.5)))), + } + + def resolve_number_and_symbol( + self, quantum_number: str, type: str, logger: BoundLogger + ) -> None: + """ + Resolves the quantum number from the `self._orbitals_map` on the passed type. `type` can be either + 'number' or 'symbol'. If the quantum type is not found, then the countertype + (e.g., type = 'number' => countertype = 'symbol') is used to resolve it. + + Args: + quantum_number (str): The quantum number to resolve. Can be 'l', 'ml' or 'ms'. + type (str): The type of the quantum number. Can be 'number' or 'symbol'. + logger (BoundLogger): The logger to log messages. + """ + if quantum_number not in ["l", "ml", "ms"]: + logger.warning( + "The quantum number is not recognized. Try 'l', 'ml' or 'ms'." + ) + return + if type not in ["number", "symbol"]: + logger.warning("The type is not recognized. Try 'number' or 'symbol'.") + return + + _countertype_map = { + "number": "symbol", + "symbol": "number", + } + # Check if quantity already exists + quantity = getattr(self, f"{quantum_number}_quantum_{type}") + if quantity or quantity is None: + return + + # If not, check whether the countertype exists + other_quantity = getattr( + self, f"{quantum_number}_quantum_{_countertype_map[type]}" + ) + if not other_quantity: + logger.warning( + f"Could not find the {quantum_number}_quantum_{type} countertype {_countertype_map[type]}." + ) + return + + # If the counterpart exists, then resolve the quantity from the orbitals_map + quantity = self._orbitals_map.get(f"{quantum_number}_{type}s", {}).get( + other_quantity + ) + setattr(self, f"{quantum_number}_quantum_{type}", quantity) + + def resolve_degeneracy(self) -> Optional[int]: + """ + Resolves the degeneracy of the orbital state. If `j_quantum_number` is not defined, then the + degeneracy is computed from the `l_quantum_number` and `ml_quantum_number`. If `j_quantum_number` + is defined, then the degeneracy is computed from the `j_quantum_number` and `mj_quantum_number`. + There is a factor of 2 included if `ms_quantum_number` is not defined (for orbitals which + are spin-degenerated). + + Returns: + (Optional[int]): The degeneracy of the orbital state. + """ + if self.l_quantum_number and self.ml_quantum_number is None: + if self.ms_quantum_number: + degeneracy = 2 * self.l_quantum_number + 1 + else: + degeneracy = 2 * (2 * self.l_quantum_number + 1) + elif self.l_quantum_number and self.ml_quantum_number: + if self.ms_quantum_number: + degeneracy = 1 + else: + degeneracy = 2 + elif self.j_quantum_number is not None: + degeneracy = 0 + for jj in self.j_quantum_number: + if self.mj_quantum_number is not None: + mjs = RussellSaundersState.generate_MJs( + self.j_quantum_number[0], rising=True + ) + degeneracy += len( + [mj for mj in mjs if mj in self.mj_quantum_number] + ) + else: + degeneracy += RussellSaundersState(jj, 1).degeneracy + return degeneracy + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Resolving the quantum numbers and symbols if not available + for quantum_number in ["l", "ml", "ms"]: + for type in ["number", "symbol"]: + self.resolve_number_and_symbol(quantum_number, type, logger) + + # Resolve the degeneracy + self.degeneracy = self.resolve_degeneracy() + + +class CoreHole(ArchiveSection): + """ + A base section used to define the core-hole state of an atom by referencing the `OrbitalsState` + section where the core-hole was generated. + """ + + orbital_ref = Quantity( + type=OrbitalsState, + description=""" + Reference to the OrbitalsState section that is used as a basis to obtain the `CoreHole` section. + """, + a_eln=ELNAnnotation(component="ReferenceEditQuantity"), + ) + + n_excited_electrons = Quantity( + type=np.float64, + description=""" + The electron charge excited for modelling purposes. This is a number between 0 and 1 (Janak state). + If `dscf_state` is set to 'initial', then this quantity is set to None (but assumed to be excited + to an excited state). + """, + ) + + dscf_state = Quantity( + type=MEnum("initial", "final"), + description=""" + Tag used to identify the role in the workflow of the same name. Allowed values are 'initial' + (not to be confused with the _initial-state approximation_) and 'final'. If 'initial' + is used, then `n_excited_electrons` is set to None and the `orbital_ref.degeneracy` is + set to 1. + """, + ) + + def resolve_occupation(self, logger: BoundLogger) -> None: + """ + Resolves the occupation of the orbital state. The occupation is resolved from the degeneracy + and the number of excited electrons. + + Args: + logger (BoundLogger): The logger to log messages. + """ + if self.orbital_ref is None or self.n_excited_electrons is None: + logger.warning( + "Cannot resolve occupation without `orbital_ref` or `n_excited_electrons`." + ) + return + if self.orbital_ref.occupation is None: + degeneracy = self.orbital_ref.resolve_degeneracy() + if degeneracy is None: + logger.warning("Cannot resolve occupation without `degeneracy`.") + return + self.orbital_ref.occupation = degeneracy - self.n_excited_electrons + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Check if n_excited_electrons is between 0 and 1 + if 0.0 <= self.n_excited_electrons <= 1.0: + logger.error("Number of excited electrons must be between 0 and 1.") + + # If dscf_state is 'initial', then n_excited_electrons is set to 0 + if self.dscf_state == "initial": + self.n_excited_electrons = None + self.degeneracy = 1 + + # Resolve the occupation of the active orbital state + if self.orbital_ref is not None and self.n_excited_electrons: + if self.orbital_ref.occupation is None: + self.resolve_occupation(logger) + + +class HubbardInteractions(ArchiveSection): + """ + A base section to define the Hubbard interactions of the system. + """ + + n_orbitals = Quantity( + type=np.int32, + description=""" + Number of orbitals used to define the Hubbard interactions. + """, + ) + + orbitals_ref = Quantity( + type=OrbitalsState, + shape=["n_orbitals"], + description=""" + Reference to the `OrbitalsState` sections that are used as a basis to obtain the Hubbard + interaction matrices. + """, + ) + + u_matrix = Quantity( + type=np.float64, + shape=["n_orbitals", "n_orbitals"], + unit="joule", + description=""" + Value of the local Hubbard interaction matrix. The order of the rows and columns coincide + with the elements in `orbital_ref`. + """, + ) + + u_interaction = Quantity( + type=np.float64, + unit="joule", + description=""" + Value of the (intraorbital) Hubbard interaction + """, + a_eln=ELNAnnotation(component="NumberEditQuantity"), + ) + + j_hunds_coupling = Quantity( + type=np.float64, + unit="joule", + description=""" + Value of the (interorbital) Hund's coupling. + """, + a_eln=ELNAnnotation(component="NumberEditQuantity"), + ) + + u_interorbital_interaction = Quantity( + type=np.float64, + unit="joule", + description=""" + Value of the (interorbital) Coulomb interaction. In rotational invariant systems, + u_interorbital_interaction = u_interaction - 2 * j_hunds_coupling. + """, + a_eln=ELNAnnotation(component="NumberEditQuantity"), + ) + + j_local_exchange_interaction = Quantity( + type=np.float64, + unit="joule", + description=""" + Value of the exchange interaction. In rotational invariant systems, j_local_exchange_interaction = j_hunds_coupling. + """, + a_eln=ELNAnnotation(component="NumberEditQuantity"), + ) + + u_effective = Quantity( + type=np.float64, + unit="joule", + description=""" + Value of the effective U parameter (u_interaction - j_local_exchange_interaction). + """, + ) + + slater_integrals = Quantity( + type=np.float64, + shape=[3], + unit="joule", + description=""" + Value of the Slater integrals [F0, F2, F4] in spherical harmonics used to derive + the local Hubbard interactions: + + u_interaction = ((2.0 / 7.0) ** 2) * (F0 + 5.0 * F2 + 9.0 * F4) / (4.0*np.pi) + + u_interorbital_interaction = ((2.0 / 7.0) ** 2) * (F0 - 5.0 * F2 + 3.0 * 0.5 * F4) / (4.0*np.pi) + + j_hunds_coupling = ((2.0 / 7.0) ** 2) * (5.0 * F2 + 15.0 * 0.25 * F4) / (4.0*np.pi) + + See e.g., Elbio Dagotto, Nanoscale Phase Separation and Colossal Magnetoresistance, + Chapter 4, Springer Berlin (2003). + """, + ) + + double_counting_correction = Quantity( + type=str, + description=""" + Name of the double counting correction algorithm applied. + """, + a_eln=ELNAnnotation(component="StringEditQuantity"), + ) + + def resolve_u_interactions(self, logger: BoundLogger) -> Optional[tuple]: + """ + Resolves the Hubbard interactions (u_interaction, u_interorbital_interaction, j_hunds_coupling) + from the Slater integrals (F0, F2, F4). + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[tuple]): The Hubbard interactions (u_interaction, u_interorbital_interaction, j_hunds_coupling). + """ + if self.slater_integrals is None or len(self.slater_integrals) == 3: + logger.warning( + "Could not find `slater_integrals` or the length is not three." + ) + return None + f0 = self.slater_integrals[0] + f2 = self.slater_integrals[1] + f4 = self.slater_integrals[2] + u_interaction = ( + ((2.0 / 7.0) ** 2) + * (f0 + 5.0 * f2 + 9.0 * f4) + / (4.0 * np.pi) + * ureg("joule") + ) + u_interorbital_interaction = ( + ((2.0 / 7.0) ** 2) + * (f0 - 5.0 * f2 + 3.0 * f4 / 2.0) + / (4.0 * np.pi) + * ureg("joule") + ) + j_hunds_coupling = ( + ((2.0 / 7.0) ** 2) + * (5.0 * f2 + 15.0 * f4 / 4.0) + / (4.0 * np.pi) + * ureg("joule") + ) + return u_interaction, u_interorbital_interaction, j_hunds_coupling + + def resolve_u_effective(self, logger: BoundLogger) -> Optional[np.float64]: + """ + Resolves the effective U parameter (u_interaction - j_local_exchange_interaction). + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[np.float64]): The effective U parameter. + """ + if self.u_interaction is None or self.j_local_exchange_interaction is None: + logger.warning( + "Could not find `HubbardInteractions.u_interaction` or `HubbardInteractions.j_local_exchange_interaction`." + ) + return None + return self.u_interaction - self.j_local_exchange_interaction + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Obtain (u, up, j_hunds_coupling) from slater_integrals + if ( + self.u_interaction is None + and self.u_interorbital_interaction is None + and self.j_hunds_coupling is None + ): + ( + self.u_interaction, + self.u_interorbital_interaction, + self.j_hunds_coupling, + ) = self.resolve_u_interactions(logger) + + # If u_effective is not available, calculate it + if self.u_effective is None: + self.u_effective = self.resolve_u_effective(logger) + + # Check if length of `orbitals_ref` is the same as the length of `umn`: + if self.u_matrix is not None and self.orbitals_ref is not None: + if len(self.u_matrix) != len(self.orbitals_ref): + logger.error( + "The length of `HubbardInteractions.u_matrix` does not coincide with length of `HubbardInteractions.orbitals_ref`." + ) + + +class AtomsState(ArchiveSection): + """ + A base section to define each atom state information. + """ + + # TODO check what happens with ghost atoms that can have `chemical_symbol='X'` + chemical_symbol = Quantity( + type=MEnum(ase.data.chemical_symbols[1:]), + description=""" + Symbol of the element, e.g. 'H', 'Pb'. This quantity is equivalent to `atomic_numbers`. + """, + ) + + atomic_number = Quantity( + type=np.int32, + description=""" + Atomic number Z. This quantity is equivalent to `chemical_symbol`. + """, + ) + + orbitals_state = SubSection(sub_section=OrbitalsState.m_def, repeats=True) + + charge = Quantity( + type=np.int32, + default=0, + description=""" + Charge of the atom. It is defined as the number of extra electrons or holes in the + atom. If the atom is neutral, charge = 0 and the summation of all (if available) the`OrbitalsState.occupation` + coincides with the `atomic_number`. Otherwise, charge can be any positive integer (+1, +2...) + for cations or any negative integer (-1, -2...) for anions. + + Note: for `CoreHole` systems we do not consider the charge of the atom even if + we do not store the final `OrbitalsState` where the electron was excited to. + """, + a_eln=ELNAnnotation(component="NumberEditQuantity"), + ) + + core_hole = SubSection(sub_section=CoreHole.m_def, repeats=False) + + hubbard_interactions = SubSection( + sub_section=HubbardInteractions.m_def, repeats=False + ) + + def resolve_chemical_symbol_and_number(self, logger: BoundLogger) -> None: + """ + Resolves the chemical symbol from the atomic number and viceversa. + + Args: + logger (BoundLogger): The logger to log messages. + """ + f = lambda x: tuple(map(bool, x)) + if f((self.chemical_symbol, self.atomic_number)) == f((None, not None)): + try: + self.chemical_symbol = ase.data.chemical_symbols[self.atomic_number] + except IndexError: + logger.error( + "The `AtomsState.atomic_number` is out of range of the periodic table." + ) + elif f((self.chemical_symbol, self.atomic_number)) == f((not None, None)): + try: + self.atomic_number = ase.data.atomic_numbers[self.chemical_symbol] + except IndexError: + logger.error( + "The `AtomsState.chemical_symbol` is not recognized in the periodic table." + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + # Get chemical_symbol from atomic_number and viceversa + self.resolve_chemical_symbol_and_number(logger) diff --git a/simulationdataschema/general.py b/simulationdataschema/general.py index ed1937a0..62f37ca8 100644 --- a/simulationdataschema/general.py +++ b/simulationdataschema/general.py @@ -193,7 +193,8 @@ def normalize(self, archive, logger) -> None: logger.error("No system information reported.") return system_ref = self.model_system[-1] - system_ref.is_representative = True + # * We define is_representative in the parser + # system_ref.is_representative = True self.m_cache["system_ref"] = system_ref # Setting up the `branch_depth` in the parent-child tree diff --git a/simulationdataschema/model_method.py b/simulationdataschema/model_method.py index 26d0bda1..30a3ca68 100644 --- a/simulationdataschema/model_method.py +++ b/simulationdataschema/model_method.py @@ -47,6 +47,6 @@ class ModelMethod(ArchiveSection): normalizer_level = 1 - def normalize(self, archive, logger): + def normalize(self, archive, logger) -> None: super().normalize(archive, logger) self.logger = logger diff --git a/simulationdataschema/model_system.py b/simulationdataschema/model_system.py index 8f3b55f9..874a0ad6 100644 --- a/simulationdataschema/model_system.py +++ b/simulationdataschema/model_system.py @@ -37,7 +37,7 @@ import re import numpy as np import ase -from typing import Union, Tuple +from typing import Tuple, Optional from structlog.stdlib import BoundLogger from matid import SymmetryAnalyzer, Classifier # pylint: disable=import-error @@ -55,19 +55,18 @@ from nomad.units import ureg from nomad.atomutils import Formula, get_normalized_wyckoff, search_aflow_prototype -from nomad.datamodel.data import ArchiveSection - from nomad.metainfo import Quantity, SubSection, SectionProxy, MEnum +from nomad.datamodel.data import ArchiveSection from nomad.datamodel.metainfo.basesections import Entity, System from nomad.datamodel.metainfo.annotations import ELNAnnotation +from .atoms_state import AtomsState from .utils import get_sibling_section class GeometricSpace(Entity): """ - A base section to define geometrical space-related entities, as well as their coordinates - system of reference information. + A base section to define geometrical space-related entities. """ length_vector_a = Quantity( @@ -192,11 +191,14 @@ class GeometricSpace(Entity): """, ) - def get_geometric_space_for_atomic_cell(self, logger): + def get_geometric_space_for_atomic_cell(self, logger: BoundLogger) -> None: """ Get the real space parameters for the atomic cell using ASE. + + Args: + logger (BoundLogger): The logger to log messages. """ - atoms = self.to_ase_atoms(logger) + atoms = self.to_ase_atoms(logger) # function defined in AtomicCell cell = atoms.get_cell() self.length_vector_a, self.length_vector_b, self.length_vector_c = ( cell.lengths() * ureg.angstrom @@ -207,16 +209,19 @@ def get_geometric_space_for_atomic_cell(self, logger): self.volume = cell.volume * ureg.angstrom**3 def normalize(self, archive, logger) -> None: - # We try to get ASE atoms from the cache of AtomicCell - ase_atoms = self.to_ase_atoms(logger) - if ase_atoms is not None: + # Skip normalization for `Entity` + try: self.get_geometric_space_for_atomic_cell(logger) + except Exception: + logger.warning( + "Could not extract the geometric space information from ASE Atoms object.", + ) + return -class AtomicCell(GeometricSpace): +class Cell(GeometricSpace): """ - A base section used to specify the atomic cell quantities (labels, positions) of a system - at a given moment in time. + A base section used to specify the cell quantities of a system at a given moment in time. """ type = Quantity( @@ -229,29 +234,16 @@ class AtomicCell(GeometricSpace): """, ) - # TODO re-define `labels` and `atomic_numbers` to handle toy-models and real chemical - # elements, and also define AtomType for these and AtomParameters - labels = Quantity( - type=str, - shape=["*"], - description=""" - List containing the labels of the atomic species in the system at the different positions - of the structure. It refers to a chemical element as defined in the periodic table, - e.g., 'H', 'O', 'Pt'. This quantity is equivalent to atomic_numbers. - """, - ) - - atomic_numbers = Quantity( + n_cell_points = Quantity( type=np.int32, - shape=["*"], description=""" - List of atomic numbers Z. This quantity is equivalent to labels. + Number of cell points. """, ) positions = Quantity( type=np.float64, - shape=["*", 3], + shape=["n_cell_points", 3], unit="meter", description=""" Positions of all the atoms in Cartesian coordinates. @@ -260,13 +252,13 @@ class AtomicCell(GeometricSpace): velocities = Quantity( type=np.float64, - shape=["*", 3], + shape=["n_cell_points", 3], unit="meter / second", description=""" Velocities of the atoms. It is the change in cartesian coordinates of the atom position with time. """, - ) # ? what about forces, stress + ) lattice_vectors = Quantity( type=np.float64, @@ -309,27 +301,45 @@ class AtomicCell(GeometricSpace): """, ) + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + +class AtomicCell(Cell): + """ + A base section used to specify the atomic cell information of a system. + """ + + atoms_state = SubSection(sub_section=AtomsState.m_def, repeats=True) + + n_atoms = Quantity( + type=np.int32, + description=""" + Number of atoms in the atomic cell. + """, + ) + equivalent_atoms = Quantity( type=np.int32, - shape=["*"], + shape=["n_atoms"], description=""" - List of equivalent atoms as defined in labels. If no equivalent atoms are found, + List of equivalent atoms as defined in `atoms`. If no equivalent atoms are found, then the list is simply the index of each element, e.g.: - [0, 1, 2, 3] all four atoms are non-equivalent. - [0, 0, 0, 3] three equivalent atoms and one non-equivalent. """, ) + # ! improve description and clarify whether this belongs to `Symmetry` with @lauri-codes wyckoff_letters = Quantity( type=str, - shape=["*"], - # TODO improve description + shape=["n_atoms"], description=""" - Wyckoff letters associated with each atom position. + Wyckoff letters associated with each atom. """, ) - def to_ase_atoms(self, logger: BoundLogger) -> Union[ase.Atoms, None]: + def to_ase_atoms(self, logger: BoundLogger) -> Optional[ase.Atoms]: """ Generates an ASE Atoms object with the most basic information from the parsed `AtomicCell` section (labels, periodic_boundary_conditions, positions, and lattice_vectors). @@ -338,10 +348,11 @@ def to_ase_atoms(self, logger: BoundLogger) -> Union[ase.Atoms, None]: logger (BoundLogger): The logger to log messages. Returns: - Union[ase.Atoms, None]: The ASE Atoms object with the basic information from the `AtomicCell`. + (Optional[ase.Atoms]): The ASE Atoms object with the basic information from the `AtomicCell`. """ # Initialize ase.Atoms object with labels - ase_atoms = ase.Atoms(symbols=self.labels) + atoms_labels = [atom_state.chemical_symbol for atom_state in self.atoms_state] + ase_atoms = ase.Atoms(symbols=atoms_labels) # PBC if self.periodic_boundary_conditions is None: @@ -353,9 +364,9 @@ def to_ase_atoms(self, logger: BoundLogger) -> Union[ase.Atoms, None]: # Positions (ensure they are parsed) if self.positions is not None: - if len(self.positions) != len(self.labels): + if len(self.positions) != len(self.atoms_state): logger.error( - "Length of `AtomicCell.positions` does not coincide with the length of the `AtomicCell.labels`." + "Length of `AtomicCell.positions` does not coincide with the length of the `AtomicCell.atoms_state`." ) return None ase_atoms.set_positions(self.positions.to("angstrom").magnitude) @@ -375,45 +386,16 @@ def to_ase_atoms(self, logger: BoundLogger) -> Union[ase.Atoms, None]: return ase_atoms - def normalize(self, archive, logger): - # If the labels and atomic_numbers are not specified, return with an error - if self.labels is None and self.atomic_numbers is None: - logger.error( - "Could not read `AtomicCell.labels` or `AtomicCell.atomic_numbers`." - ) - return - atomic_labels = self.labels - atomic_numbers = self.atomic_numbers - - # Labels - if atomic_labels is None and atomic_numbers is not None: - try: - atomic_labels = [ - ase.data.chemical_symbols[number] for number in atomic_numbers - ] - except IndexError: - logger.error( - "The `AtomicCell.atomic_numbers` are out of range of the periodic table." - ) - return - self.labels = atomic_labels - - # Use ASE Atoms functionalities to extract information about the AtomicCell - ase_atoms = self.to_ase_atoms(logger) - - # Atomic numbers - if atomic_labels is not None and atomic_numbers is None: - atomic_numbers = ase_atoms.get_atomic_numbers() - self.atomic_numbers = atomic_numbers - - # We then normalize `GeometricSpace` + def normalize(self, archive, logger) -> None: super().normalize(archive, logger) class Symmetry(ArchiveSection): """ - A base section used to specify the symmetry of the `AtomicCell`. This information can - be extracted via normalization using the MatID package, if `AtomicCell` is specified. + A base section used to specify the symmetry of the `AtomicCell`. + + Note: this information can be extracted via normalization using the MatID package, if `AtomicCell` + is specified. """ bravais_lattice = Quantity( @@ -519,7 +501,7 @@ class Symmetry(ArchiveSection): def resolve_analyzed_atomic_cell( self, symmetry_analyzer: SymmetryAnalyzer, cell_type: str, logger: BoundLogger - ) -> Union[AtomicCell, None]: + ) -> Optional[AtomicCell]: """ Resolves the `AtomicCell` section from the `SymmetryAnalyzer` object and the cell_type (primitive or conventional). @@ -530,7 +512,7 @@ def resolve_analyzed_atomic_cell( logger (BoundLogger): The logger to log messages. Returns: - Union[AtomicCell, None]: The resolved `AtomicCell` section or None if the cell_type + (Optional[AtomicCell]): The resolved `AtomicCell` section or None if the cell_type is not recognized. """ if cell_type not in ["primitive", "conventional"]: @@ -551,11 +533,12 @@ def resolve_analyzed_atomic_cell( atomic_cell = AtomicCell(type=cell_type) atomic_cell.lattice_vectors = cell * ureg.angstrom + for label, atomic_number in zip(labels, atomic_numbers): + atom_state = AtomsState(chemical_symbol=label, atomic_number=atomic_number) + atomic_cell.m_add_sub_section(AtomicCell.atoms_state, atom_state) atomic_cell.positions = ( positions * ureg.angstrom ) # ? why do we need to pass units - atomic_cell.labels = labels - atomic_cell.atomic_numbers = atomic_numbers atomic_cell.wyckoff_letters = wyckoff atomic_cell.equivalent_atoms = equivalent_atoms atomic_cell.get_geometric_space_for_atomic_cell(logger) @@ -563,7 +546,7 @@ def resolve_analyzed_atomic_cell( def resolve_bulk_symmetry( self, original_atomic_cell: AtomicCell, logger: BoundLogger - ) -> Tuple[Union[AtomicCell, None], Union[AtomicCell, None]]: + ) -> Tuple[Optional[AtomicCell], Optional[AtomicCell]]: """ Resolves the symmetry of the material being simulated using MatID and the originally parsed data under original_atomic_cell. It generates two other @@ -575,8 +558,8 @@ def resolve_bulk_symmetry( uses to in MatID.SymmetryAnalyzer(). logger (BoundLogger): The logger to log messages. Returns: - primitive_atomic_cell (AtomicCell): The primitive `AtomicCell` section. - conventional_atomic_cell (AtomicCell): The standarized `AtomicCell` section. + primitive_atomic_cell (Optional[AtomicCell]): The primitive `AtomicCell` section. + conventional_atomic_cell (Optional[AtomicCell]): The standarized `AtomicCell` section. """ symmetry = {} try: @@ -612,21 +595,27 @@ def resolve_bulk_symmetry( original_atomic_cell.wyckoff_letters = original_wyckoff original_atomic_cell.equivalent_atoms = original_equivalent_atoms - # Populating the primitive AtomicCell information + # Populating the primitive AtomState information primitive_atomic_cell = self.resolve_analyzed_atomic_cell( symmetry_analyzer, "primitive", logger ) - # Populating the conventional Atoms information + # Populating the conventional AtomState information conventional_atomic_cell = self.resolve_analyzed_atomic_cell( symmetry_analyzer, "conventional", logger ) # Getting prototype_formula, prototype_aflow_id, and strukturbericht designation from # standarized Wyckoff numbers and the space group number - if symmetry.get("space_group_number"): + if ( + symmetry.get("space_group_number") + and conventional_atomic_cell.atoms_state is not None + ): # Resolve atomic_numbers and wyckoff letters from the conventional cell - conventional_num = conventional_atomic_cell.atomic_numbers + conventional_num = [ + atom_state.atomic_number + for atom_state in conventional_atomic_cell.atoms_state + ] conventional_wyckoff = conventional_atomic_cell.wyckoff_letters # Normalize wyckoff letters norm_wyckoff = get_normalized_wyckoff( @@ -636,8 +625,11 @@ def resolve_bulk_symmetry( symmetry.get("space_group_number"), norm_wyckoff ) strukturbericht = aflow_prototype.get("Strukturbericht Designation") - if strukturbericht != "None": - strukturbericht = re.sub("[$_{}]", "", strukturbericht) + strukturbericht = ( + re.sub("[$_{}]", "", strukturbericht) + if strukturbericht != "None" + else None + ) prototype_aflow_id = aflow_prototype.get("aflow_prototype_id") prototype_formula = aflow_prototype.get("Prototype") # Adding these to the symmetry dictionary for later assignement @@ -651,7 +643,7 @@ def resolve_bulk_symmetry( return primitive_atomic_cell, conventional_atomic_cell - def normalize(self, archive, logger): + def normalize(self, archive, logger) -> None: atomic_cell = get_sibling_section( section=self, sibling_section_name="atomic_cell", logger=logger ) @@ -727,7 +719,7 @@ class ChemicalFormula(ArchiveSection): """, ) - def resolve_chemical_formulas(self, formula: Formula): + def resolve_chemical_formulas(self, formula: Formula) -> None: """ Resolves the chemical formulas of the `ModelSystem` in different formats. @@ -740,7 +732,7 @@ def resolve_chemical_formulas(self, formula: Formula): self.hill = formula.format("hill") self.anonymous = formula.format("anonymous") - def normalize(self, archive, logger): + def normalize(self, archive, logger) -> None: atomic_cell = get_sibling_section( section=self, sibling_section_name="atomic_cell", logger=logger ) @@ -761,37 +753,44 @@ def normalize(self, archive, logger): class ModelSystem(System): + # ! Work in this description with new `AtomState` """ - Model system used as an input for the computation. It inherits from `System` where a set - of sub-sections for the `elemental_composition` is defined. + Model system used as an input for the computation. - We also defined: + We define: - `name` refers to all the verbose and user-dependent naming in ModelSystem, - `type` refers to the type of the ModelSystem (atom, bulk, surface, etc.), - - `dimensionality` refers to the dimensionality of the ModelSystem (0D, 1D, 2D, 3D), + - `dimensionality` refers to the dimensionality of the ModelSystem (0, 1, 2, 3), - If the ModelSystem `is_representative`, the normalization occurs. The time evolution of - the system is encoded on the fact that ModelSystem is a list under Simulation, and for - each element of that list, `time_step` can be defined. + If the ModelSystem `is_representative`, proceeds with normalization. The time evolution of the + ModelSystem is stored in a `list` format under `Simulation`, and for each element of that list, + `time_step` can be defined. It is composed of the sub-sections: - `AtomicCell` containing the information of the atomic structure, - - `Symmetry` containing the information of the (standarized) atomic cell symmetry + - `Symmetry` containing the information of the (conventional) atomic cell symmetry in bulk ModelSystem, - `ChemicalFormula` containing the information of the chemical formulas in different formats. - This class nest over itself (with the proxy in `model_system`) to define different + This class nest over itself (with the section proxy in `model_system`) to define different parent-child system trees. The quantities `branch_label`, `branch_depth`, `atom_indices`, and `bond_list` are used to define the parent-child tree. The normalization is ran in the following order: - 1. `AtomicCell.normalize()` from `atomic_cell`, - 2. `ModelSystem.normalize()` in this class, - 3. `Symmetry.normalize()` is called within this class normalization, - 4. `ChemicalFormula.normalize()` is called within this class normalization. + 1. `OrbitalsState.normalize()` in atoms_state.py under `AtomsState` + 2. `CoreHoleState.normalize()` in atoms_state.py under `AtomsState` + 3. `HubbardInteractions.normalize()` in atoms_state.py under `AtomsState` + 4. `AtomsState.normalize()` in atoms_state.py + 5. `AtomicCell.normalize()` in atomic_cell.py + 6. `Symmetry.normalize()` in this class + 7. `ChemicalFormula.normalize()` in this class + 8. `ModelSystem.normalize()` in this class - Examples: + Note: `normalize()` can be called at any time for each of the classes without being re-triggered + by the NOMAD normalization. + + Examples for the parent-child hierarchical trees: - Example 1, a crystal Si has: 3 AtomicCell sections (named 'original', 'primitive', and 'conventional'), 1 Symmetry section, and 0 nested ModelSystem trees. @@ -866,6 +865,7 @@ class ModelSystem(System): """, ) + # ? Check later when implementing `Outputs` if this quantity needs to be extended time_step = Quantity( type=np.int32, description=""" @@ -883,7 +883,6 @@ class ModelSystem(System): branch_label = Quantity( type=str, - shape=[], description=""" Label of the specific branch in the hierarchical `ModelSystem` tree. """, @@ -896,7 +895,6 @@ class ModelSystem(System): """, ) - # TODO add method to resolve labels and positions from the parent AtomicCell atom_indices = Quantity( type=np.int32, shape=["*"], @@ -909,9 +907,9 @@ class ModelSystem(System): """, ) + # TODO improve description and add an example using the case in atom_indices bond_list = Quantity( type=np.int32, - # TODO improve description and add an example using the case in atom_indices description=""" List of pairs of atom indices corresponding to bonds (e.g., as defined by a force field) within this atoms_group. @@ -921,10 +919,10 @@ class ModelSystem(System): model_system = SubSection(sub_section=SectionProxy("ModelSystem"), repeats=True) def resolve_system_type_and_dimensionality( - self, ase_atoms: ase.Atoms + self, ase_atoms: ase.Atoms, logger: BoundLogger ) -> Tuple[str, int]: """ - Determine the ModelSystem.type and ModelSystem.dimensionality using MatID classification analyzer: + Resolves the `ModelSystem.type` and `ModelSystem.dimensionality` using `MatID` classification analyzer: - https://singroup.github.io/matid/tutorials/classification.html @@ -947,7 +945,7 @@ def resolve_system_type_and_dimensionality( ) classification = classifier.classify(ase_atoms) except Exception as e: - self.logger.warning( + logger.warning( "MatID system classification failed.", exc_info=e, error=str(e) ) return system_type, dimensionality @@ -971,15 +969,14 @@ def resolve_system_type_and_dimensionality( system_type = "2D" dimensionality = 2 else: - self.logger.info( + logger.info( "ModelSystem.type and dimensionality analysis not run due to large system size." ) return system_type, dimensionality - def normalize(self, archive, logger): + def normalize(self, archive, logger) -> None: super().normalize(archive, logger) - self.logger = logger # We don't need to normalize if the system is not representative if not self.is_representative: @@ -988,8 +985,7 @@ def normalize(self, archive, logger): # Extracting ASE Atoms object from the originally parsed AtomicCell section if self.atomic_cell is None: self.logger.warning( - "Could not find the originally parsed atomic system. " - "Symmetry and ChemicalFormula extraction is thus not run." + "Could not find the originally parsed atomic system. `Symmetry` and `ChemicalFormula` extraction is thus not run." ) return self.atomic_cell[0].type = "original" @@ -1005,7 +1001,7 @@ def normalize(self, archive, logger): ( self.type, self.dimensionality, - ) = self.resolve_system_type_and_dimensionality(ase_atoms) + ) = self.resolve_system_type_and_dimensionality(ase_atoms, logger) # Creating and normalizing Symmetry section if self.type == "bulk" and self.symmetry is not None: sec_symmetry = self.m_create(Symmetry) diff --git a/simulationdataschema/outputs.py b/simulationdataschema/outputs.py index f0230b9f..6b1dfbd1 100644 --- a/simulationdataschema/outputs.py +++ b/simulationdataschema/outputs.py @@ -47,6 +47,6 @@ class Outputs(ArchiveSection): normalizer_level = 2 - def normalize(self, archive, logger): + def normalize(self, archive, logger) -> None: super().normalize(archive, logger) self.logger = logger diff --git a/simulationdataschema/utils/__init__.py b/simulationdataschema/utils/__init__.py index 0cd04b48..d41b8bee 100644 --- a/simulationdataschema/utils/__init__.py +++ b/simulationdataschema/utils/__init__.py @@ -16,4 +16,4 @@ # See the License for the specific language governing permissions and # limitations under the License. -from .utils import get_sibling_section +from .utils import get_sibling_section, RussellSaundersState diff --git a/simulationdataschema/utils/utils.py b/simulationdataschema/utils/utils.py index 398d9aaf..e642848c 100644 --- a/simulationdataschema/utils/utils.py +++ b/simulationdataschema/utils/utils.py @@ -17,6 +17,7 @@ # limitations under the License. # +from math import factorial from typing import Optional from structlog.stdlib import BoundLogger @@ -64,3 +65,47 @@ def get_sibling_section( return None return sibling_section[index_sibling] return sibling_section + + +# ? Check if this utils deserves its own file after extending it +class RussellSaundersState: + @classmethod + def generate_Js(cls, J1: float, J2: float, rising=True): + J_min, J_max = sorted([abs(J1), abs(J2)]) + generator = range( + int(J_max - J_min) + 1 + ) # works for both for fermions and bosons + if rising: + for jj in generator: + yield J_min + jj + else: + for jj in generator: + yield J_max - jj + + @classmethod + def generate_MJs(cls, J, rising=True): + generator = range(int(2 * J + 1)) + if rising: + for m in generator: + yield -J + m + else: + for m in generator: + yield J - m + + def __init__(self, *args, **kwargs): + self.J = kwargs.get("J") + if self.J is None: + raise TypeError + self.occupation = kwargs.get("occ") + if self.occupation is None: + raise TypeError + + @property + def multiplicity(self): + return 2 * self.J + 1 + + @property + def degeneracy(self): + return factorial(self.multiplicity) / ( + factorial(self.multiplicity - self.occupation) * factorial(self.occupation) + )