H5md refactor (#196)

* add velocities and particle_state labels

* working h5md basic system

* added some ModelSystem functions for molecules

* add some system types

* fixed tests

* cleaned comments

* returned current docstrings that diverged somehow

* ruff

* refined is_molecule, updated tests

* remove duplicate is_molecule

* fix editing mistakes

* switched to get_symbols() but left the old function commented for now

* re-uncomment get_chemical_symbols() for now

Author: JFRudzinski

src/nomad_simulations/schema_packages/atoms_state.py

@@ -551,7 +551,12 @@ class ParticleState(Entity):
     This can be extended to include any common quantities in the future.
     """
 
-    pass
+    label = Quantity(
+        type=str,
+        description="""
+        User- or program-package-defined identifier for this particle.
+        """,
+    )
 
 
 class AtomsState(ParticleState):
@@ -660,5 +665,83 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
             self.atomic_number = self.resolve_atomic_number(logger=logger)
 
 
-class CoarseGrainedState(ParticleState):
-    pass
+class CGBeadState(ParticleState):
+    """
+    A section to define coarse-grained bead state information.
+    """
+
+    # ? What do we want to qualify as type identifier? What safety checks do we need?
+    bead_symbol = Quantity(
+        type=str,
+        description="""
+        Symbol(s) describing the (base) CG particle type. Equivalent to chemical_symbol
+        for atomic elements.
+        """,
+    )
+
+    label = Quantity(
+        type=str,
+        description="""
+        User- or program-package-defined identifier for this bead site.
+        This could be used to store primary FF labels in cases where only a
+        secondary specification is required. Otherwise, `alt_labels` are
+        used to document more complex bead identifiers, e.g., bead interactions based
+        on connectivity.
+        """,
+    )
+
+    alt_labels = Quantity(
+        type=str,
+        shape=['*'],
+        description="""
+        A list of bead labels for multifaceted bead characterization.
+        """,
+    )
+
+    mass = Quantity(
+        type=np.float64,
+        unit='kg',
+        description="""
+        Total mass of the particle.
+        """,
+    )
+
+    charge = Quantity(
+        type=np.float64,
+        unit='coulomb',
+        description="""
+        Total charge of the particle.
+        """,
+    )
+
+    # Other possible quantities
+    #     diameter: float
+    #         The diameter of each particle.
+    #         Default: 1.0
+    #     body: int
+    #         The composite body associated with each particle. The value -1
+    #         indicates no body.
+    #         Default: -1
+    #     moment_inertia: float
+    #         The moment_inertia of each particle (I_xx, I_yy, I_zz).
+    #         This inertia tensor is diagonal in the body frame of the particle.
+    #         The default value is for point particles.
+    #         Default: 0, 0, 0
+    #     scaled_positions: list of scaled-positions #! for cell if relevant
+    #         Like positions, but given in units of the unit cell.
+    #         Can not be set at the same time as positions.
+    #         Default: 0, 0, 0
+    #     orientation: float
+    #         The orientation of each particle. In scalar + vector notation,
+    #         this is (r, a_x, a_y, a_z), where the quaternion is q = r + a_xi + a_yj + a_zk.
+    #         A unit quaternion has the property: sqrt(r^2 + a_x^2 + a_y^2 + a_z^2) = 1.
+    #         Default: 0, 0, 0, 0
+    #     angmom: float #? for cell or here?
+    #         The angular momentum of each particle as a quaternion.
+    #         Default: 0, 0, 0, 0
+    #     image: int #! advance PBC stuff would go in cell I guess
+    #         The number of times each particle has wrapped around the box (i_x, i_y, i_z).
+    #         Default: 0, 0, 0
+
+    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
+        super().normalize(archive, logger)

src/nomad_simulations/schema_packages/model_system.py

@@ -23,6 +23,7 @@
 
 import ase
 import numpy as np
+from ase.symbols import symbols2numbers
 from matid import Classifier, SymmetryAnalyzer  # pylint: disable=import-error
 from matid.classification.classifications import (
     Atom,
@@ -51,6 +52,7 @@
 
 from nomad_simulations.schema_packages.atoms_state import (
     AtomsState,
+    CGBeadState,
     ParticleState,
 )
 from nomad_simulations.schema_packages.utils import (
@@ -844,7 +846,9 @@ class ModelSystem(System):
         type=MEnum(
             'atom',
             'active_atom',
-            'molecule / cluster',
+            'molecule',
+            'monomer',
+            'cluster',
             '1D',
             'surface',
             '2D',
@@ -911,7 +915,7 @@ class ModelSystem(System):
         type=np.int32,
         shape=['*'],
         description="""
-        Global indices of the particles that belong to this subsystem, 
+        Global indices of the particles that belong to this subsystem,
         counted from the representative (top-level) ModelSystem.
 
         **Example (SrTiO_3 primitive cell)**
@@ -938,6 +942,16 @@ class ModelSystem(System):
         """,
     )
 
+    velocities = Quantity(
+        type=np.float64,
+        shape=['*', 3],
+        unit='meter / second',
+        description="""
+            Velocities of the particles: I.e., the change in cartesian coordinates of the
+        particle position with time.
+        """,
+    )
+
     # TODO improve description and add an example
     bond_list = Quantity(
         type=np.int32,
@@ -993,7 +1007,7 @@ class ModelSystem(System):
         section_def=ParticleState.m_def,
         repeats=True,
         description="""
-        Particle state of each of the particles conforming the ModelSystem. 
+        Particle state of each of the particles conforming the ModelSystem.
         This is a list of `n_particles` elements and the order matches that of `positions`.
 
             Example
@@ -1010,6 +1024,7 @@ class ModelSystem(System):
 
     sub_systems = SubSection(sub_section=SectionProxy('ModelSystem'), repeats=True)
 
+    # TODO Will remove this after developing CGBeadState functionality further
     def get_chemical_symbols(self, logger: 'BoundLogger') -> list[str]:
         """
         Gets the chemical symbols from the particle_states that are AtomsState instances.
@@ -1028,6 +1043,55 @@ def get_chemical_symbols(self, logger: 'BoundLogger') -> list[str]:
                 chemical_symbols.append(particle_state.chemical_symbol)
         return chemical_symbols
 
+    # TODO: symbols should be a property right?
+    # ? To replace get_chemical_symbols
+    def get_symbols(self, logger: 'BoundLogger') -> list[str]:
+        """
+        Gets the symbols from the particle_states.
+        Args:
+            logger (BoundLogger): The logger to log messages.
+        Returns:
+            list: The list of symbols of the particles.
+        """
+        symbols = []
+        for particle_state in self.particle_states:
+            symbol = None
+            if isinstance(particle_state, AtomsState):
+                symbol = particle_state.chemical_symbol
+            elif isinstance(particle_state, CGBeadState):
+                symbol = particle_state.bead_symbol
+            if not symbol:
+                logger.warning('missing symbol in ParticleState.')
+            symbols.append(symbol)
+        return symbols
+
+    def are_valid_chemical_symbols(self, logger: 'BoundLogger') -> bool:
+        """
+        Validate that ASE can map all element symbols in the particle_states
+        to atomic numbers.
+        Args:
+            logger (BoundLogger): The logger to log messages.
+        Returns:
+            bool: True if all chemical symbols are valid, False otherwise.
+        """
+        symbols = self.get_symbols(logger)
+        if not symbols:
+            return False
+
+        try:
+            symbols2numbers(symbols)
+        except KeyError as e:
+            logger.error(f'Invalid chemical symbol found: {e}')
+            return False
+        return True
+
+    # atom_labels = self.traj_parser.get_atom_labels(n)
+    # if atom_labels is not None:
+    #     try:
+    #         symbols2numbers(atom_labels)
+    #     except KeyError:
+    #         atom_labels = ['X'] * len(atom_labels)
+
     def to_ase_atoms(self, logger: 'BoundLogger') -> 'Optional[ase.Atoms]':
         """
         Generates an ASE Atoms object from ModelSystem data.
@@ -1036,7 +1100,7 @@ def to_ase_atoms(self, logger: 'BoundLogger') -> 'Optional[ase.Atoms]':
           - positions from the top-level positions quantity,
           - periodic boundary conditions and lattice vectors from the first cell.
         """
-        symbols = self.get_chemical_symbols(logger)
+        symbols = self.get_symbols(logger)
         if not symbols:
             logger.error('Cannot generate ASE Atoms without chemical symbols.')
             return None
@@ -1207,3 +1271,100 @@ def is_ge_structure(self, other: 'ModelSystem') -> bool:
 
     def is_ne_structure(self, other: 'ModelSystem') -> bool:
         return not self.is_equal_structure(other)
+
+    # functions for traversing the ModelSystem hierarchy
+    def get_root_system(self) -> 'ModelSystem':
+        """
+        Traverses up the hierarchy to find the root ModelSystem.
+
+        Returns:
+            ModelSystem: The top-level (root) ModelSystem.
+        """
+        system = self
+        while isinstance(system.m_parent, ModelSystem):
+            system = system.m_parent
+        return system
+
+    # functions for working with molecules
+    def get_bond_list(self, set_local: bool = False) -> np.ndarray:
+        """
+        Retrieves the bond list for this subsystem by filtering the root bond_list
+        using the subsystem's `particle_indices`. The bond indices remain in root-level
+        coordinates (no reindexing).
+
+        Args:
+            set_local (bool): If True, sets `self.bond_list` to the filtered bonds.
+
+        Returns:
+            np.ndarray: Filtered bond list for this subsystem (root-level indices).
+        """
+
+        if not isinstance(self.m_parent, ModelSystem):  # this is the root system
+            return self.bond_list
+
+        if self.particle_indices is None:
+            return np.array([])
+
+        root = self.get_root_system()
+        if root.bond_list is None:
+            return np.array([])
+
+        indices_set = set(self.particle_indices.tolist())
+        bond_list = np.array(
+            [
+                (i, j)
+                for i, j in root.bond_list
+                if i in indices_set and j in indices_set
+            ],
+            dtype=np.int32,
+        )
+
+        if set_local:
+            self.bond_list = bond_list
+
+        return bond_list
+
+    def is_molecule(self) -> bool:
+        """
+        Checks if the current subsystem forms a contiguous and isolated molecule:
+        - All particles are connected (single connected component).
+        - No bonds connect particles inside this subsystem to particles outside it.
+
+        Returns:
+            bool: True if the subsystem is an isolated molecule, False otherwise.
+        """
+        import networkx as nx
+
+        # Internal bonds for this subsystem
+        bonds = self.get_bond_list(set_local=False)
+
+        # Handle case: no bonds
+        if bonds.size == 0:
+            return False
+
+        # Determine particle indices (fallback to range if None)
+        particle_indices = self.particle_indices
+        if particle_indices is None:
+            n_particles = (
+                len(self.positions) if self.positions is not None else self.n_particles
+            )
+            particle_indices = np.arange(n_particles, dtype=np.int32)
+
+        # --- 1. Connectivity check ---
+        graph = nx.Graph()
+        graph.add_nodes_from(particle_indices)
+        graph.add_edges_from(bonds)
+
+        if not nx.is_connected(graph):
+            return False
+
+        # --- 2. Isolation check: ensure no bonds cross subsystem boundary ---
+        root = self.get_root_system()
+        if root.bond_list is not None:
+            indices_set = set(particle_indices.tolist())
+            for i, j in root.bond_list:
+                # If exactly one endpoint is inside → cross-boundary bond
+                if (i in indices_set) ^ (j in indices_set):
+                    return False
+
+        return True