Draft: add support for nanoparticle models

CHANGELOG.md

@@ -14,6 +14,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Utility functions to load and save the new database via the pyMBE API, `pmb.save_database` and `pmb.load_database`. (#147)
 - Added functions to define particle states:  `pmb.define_particle_states`  and `pmb.define_monoprototic_particle_states`. (#147)
 - Added utility functions in `lib/handy_functions` to define residue and particle templates for aminoacids en peptides and residues: `define_protein_AA_particles`, `define_protein_AA_residues` and `define_peptide_AA_residues`. (#147)
+- Added support for nanoparticle templates and instances in the canonical storage layer via `NanoparticleTemplate` and `NanoparticleInstance`. (#148)
+- Added new API methods `pmb.define_nanoparticle` and `pmb.create_nanoparticle` to define and build nanoparticles with configurable core particles and surface site composition. (#148)
+- Added nanoparticle site-construction utilities in `pyMBE.lib.nanoparticle_tools` for spherical site distribution, patch construction, and overlap checks. (#148)
+- Added sample script `samples/nanoparticles_grxmc.py` to demonstrate nanoparticle setup and simulation workflows. (#148)
+- Added dedicated nanoparticle unit tests (`testsuite/nanoparticle_unit_tests.py`) and coverage for nanoparticle-related code paths. (#148)
 
 ## Changed
 - Create methods (`create_particle`, `create_residue`, `create_molecule`, `create_protein`, `create_hydrogel`) now raise a ValueError if no template is found for an input `name` instead than a warning. (#147)
@@ -22,10 +27,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Pka values are now stored as part of chemical reactions and no longer an attribute of particle templates. (#147)
 - Amino acid residue templates are no longer defined internally in `define_peptide` and `define_protein`. Those definitions are now exposed to the user. (#147)
 - Molecule templates now need to be defined to be used as templates for hydrogel chains in hydrogels. (#147)
+- Rigid-body setup is now integrated into nanoparticle creation, allowing nanoparticles to be initialized as rigid objects directly from `pmb.create_nanoparticle`. (#148)
+- Nanoparticle construction now supports primary/secondary site partitioning and multi-patch layouts driven by template parameters. (#148)
 
 ## Fixed
 - Wrong handling of units in `get_radius_map` when the `dimensionless` argument was triggered. (#147)
 - Utility methods `get_particle_id_map`, `calculate_HH`, `calculate_net_charge`,  `center_object_in_simulation_box` now support all template types in pyMBE, including hydrogels. Some of these methods have been renamed to expose directly in the API this change in behavior. (#147)
+- Fixed edge cases in rigid-body initialization used by nanoparticle creation to improve robustness of newly created nanoparticle objects. (#148)
 
 
 ### Removed

pyMBE/lib/handy_functions.py

@@ -342,6 +342,90 @@ def get_number_of_particles(espresso_system, ptype):
         kwargs = {"type": ptype}
     return espresso_system.number_of_particles(*args, **kwargs)
 
+def generate_lattice_positions(lattice_type, number_of_sites, lattice_constant=1.0, box_length=None, origin=None):
+    """
+    Generate lattice positions for a requested lattice type and number of sites.
+
+    Args:
+        lattice_type ('str'):
+            Lattice type identifier. Supported values are:
+            - ``"sc"``  (simple cubic)
+            - ``"bcc"`` (body-centered cubic)
+            - ``"fcc"`` (face-centered cubic)
+
+        number_of_sites ('int'):
+            Number of lattice positions to generate.
+
+        lattice_constant ('float', optional):
+            Lattice constant. Used when ``box_length`` is not provided.
+            Must be positive.
+
+        box_length ('float', optional):
+            If provided, lattice positions are fitted into a cubic box of side
+            ``box_length`` by choosing the cell spacing automatically.
+
+        origin ('list[float]', optional):
+            Origin shift applied to all generated coordinates.
+            Defaults to ``[0.0, 0.0, 0.0]``.
+
+    Returns:
+        ('list[list[float]]'):
+            List of 3D lattice positions.
+
+    Raises:
+        ValueError:
+            If ``lattice_type`` is unsupported, ``number_of_sites`` is negative,
+            or geometric inputs are invalid.
+    """
+    lattice_key = lattice_type.lower()
+    basis_map = {
+        "sc": np.array([[0.0, 0.0, 0.0]]),
+        "bcc": np.array([[0.0, 0.0, 0.0],
+                         [0.5, 0.5, 0.5]]),
+        "fcc": np.array([[0.0, 0.0, 0.0],
+                         [0.0, 0.5, 0.5],
+                         [0.5, 0.0, 0.5],
+                         [0.5, 0.5, 0.0]]),
+    }
+    if lattice_key not in basis_map:
+        raise ValueError(f"Unsupported lattice_type '{lattice_type}'. Supported values are {list(basis_map.keys())}.")
+    if number_of_sites < 0:
+        raise ValueError("number_of_sites must be a non-negative integer.")
+    if number_of_sites == 0:
+        return []
+    if origin is None:
+        origin = np.zeros(3)
+    else:
+        origin = np.array(origin, dtype=float)
+        if origin.shape != (3,):
+            raise ValueError("origin must be a 3D coordinate [x, y, z].")
+
+    points_per_cell = len(basis_map[lattice_key])
+    n_cells = int(np.ceil((number_of_sites / points_per_cell) ** (1.0 / 3.0)))
+    if n_cells <= 0:
+        n_cells = 1
+
+    if box_length is not None:
+        if box_length <= 0:
+            raise ValueError("box_length must be positive.")
+        spacing = float(box_length) / n_cells
+    else:
+        if lattice_constant <= 0:
+            raise ValueError("lattice_constant must be positive.")
+        spacing = float(lattice_constant)
+
+    basis = basis_map[lattice_key]
+    positions = []
+    for i in range(n_cells):
+        for j in range(n_cells):
+            for k in range(n_cells):
+                cell_origin = np.array([i, j, k], dtype=float) * spacing
+                for site in basis:
+                    positions.append((cell_origin + site * spacing + origin).tolist())
+                    if len(positions) == number_of_sites:
+                        return positions
+    return positions
+
 def get_residues_from_topology_dict(topology_dict, model):
     """
     Groups beads from a topology dictionary into residues and assigns residue names.
@@ -751,4 +835,4 @@ def setup_electrostatic_interactions(units, espresso_system, kT, c_salt=None, so
         espresso_system.actors.add(coulomb)
     else:
         espresso_system.electrostatics.solver = coulomb
-    logging.debug("*** Electrostatics successfully added to the system ***")
+    logging.debug("*** Electrostatics successfully added to the system ***")