diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7a1e3fad..cdbd994f 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/pyMBE/lib/handy_functions.py b/pyMBE/lib/handy_functions.py
index 33da68b0..a90c2584 100644
--- a/pyMBE/lib/handy_functions.py
+++ b/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 ***")
\ No newline at end of file
+ logging.debug("*** Electrostatics successfully added to the system ***")
diff --git a/pyMBE/lib/nanoparticle_tools.py b/pyMBE/lib/nanoparticle_tools.py
new file mode 100644
index 00000000..2e1c8d6a
--- /dev/null
+++ b/pyMBE/lib/nanoparticle_tools.py
@@ -0,0 +1,300 @@
+
+#
+# Copyright (C) 2024 pyMBE-dev team
+#
+# This file is part of pyMBE.
+#
+# pyMBE is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# pyMBE is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+import numpy as np
+import random
+from scipy.spatial import cKDTree
+
+# Distributing points evenly on the surface of a sphere
+
+def uniform_distribution_sites_on_sphere(number_of_edges=2, tolerance=1e-6):
+ """
+ Distribute points approximately uniformly on the surface of a unit sphere.
+
+ The algorithm uses iterative force-based relaxation, conceptually similar
+ to the Thomson problem.
+
+ Args:
+ number_of_edges ('int', optional):
+ Number of points to distribute on the sphere surface.
+
+ tolerance ('float', optional):
+ Convergence tolerance for the iterative relaxation.
+
+ Returns:
+ ('list[tuple[float, float, float]]'):
+ Point coordinates on the unit sphere centered at ``[0, 0, 0]``.
+
+ Notes:
+ – Thomson problem — Wikipedia (overview of the physics problem), Wikipedia.
+ – Simple schemes for uniform point distribution. Cheng Guan Koay, J Comput Sci. 2011 Dec;2(4):377–381. doi: 10.1016/j.jocs.2011.06.007.
+ """
+
+ # Generates initial configuration
+
+ random.seed(42)
+ edges = []
+ for i in range(number_of_edges):
+
+ theta = random.random() * 2*np.pi
+ phi = np.arcsin(random.random() * 2 - 1)
+ edges.append((np.cos(theta)*np.cos(phi), np.sin(theta)*np.cos(phi), np.sin(phi)))
+
+ # Iterates until fullfil the tolerance
+
+ while 1:
+ # Determine the total force acting on each point.
+ forces = []
+ for i in range(len(edges)):
+ p = edges[i]
+ f = (0,0,0)
+ ftotal = 0
+ for j in range(len(edges)):
+ if j == i: continue
+ q = edges[j]
+
+ # Find the distance vector, and its length.
+ dv = (p[0]-q[0], p[1]-q[1], p[2]-q[2])
+ dl = np.sqrt(dv[0]**2 + dv[1]**2 + dv[2]**2)
+
+ # The force vector is dv divided by dl^3. (We divide by dl once to make dv a unit vector, then by dl^2 to make its length correspond to the force.)
+ dl3 = dl ** 3
+ fv = (dv[0]/dl3, dv[1]/dl3, dv[2]/dl3)
+
+ # Add to the total force on the point p.
+ f = (f[0]+fv[0], f[1]+fv[1], f[2]+fv[2])
+
+ # Add in the forces array.
+ forces.append(f)
+
+ # Add to the running sum of the total forces/distances.
+ ftotal = ftotal + np.sqrt(f[0]**2 + f[1]**2 + f[2]**2)
+
+ # Scale the forces to ensure the points do not move too far in one go. Otherwise there will be chaotic jumping around and never any convergence.
+
+ if ftotal > 0.25:
+ fscale = 0.25 / ftotal
+ else:
+ fscale = 1
+
+ # Move each point, and normalise. While we do this, also track the distance each point ends up moving.
+
+ dist = 0
+ for i in range(len(edges)):
+ p = edges[i]
+ f = forces[i]
+ p2 = (p[0] + f[0]*fscale, p[1] + f[1]*fscale, p[2] + f[2]*fscale)
+ pl = np.sqrt(p2[0]**2 + p2[1]**2 + p2[2]**2)
+ p2 = (p2[0] / pl, p2[1] / pl, p2[2] / pl)
+ dv = (p[0]-p2[0], p[1]-p2[1], p[2]-p2[2])
+ dl = np.sqrt(dv[0]**2 + dv[1]**2 + dv[2]**2)
+ dist = dist + dl
+ edges[i] = p2
+
+ # Check for convergence and finish.
+
+ if dist < tolerance:
+ break
+
+ return edges
+
+# Auxiliary functions to calculate the patchy distribution
+
+def calculate_distance_vector_point(A,p):
+ """
+ Compute Euclidean distances between a set of 3D points and one 3D point.
+
+ Args:
+ A ('iterable'):
+ Iterable of 3D points.
+
+ p ('iterable'):
+ Reference 3D point.
+
+ Returns:
+ ('list[float]'):
+ Euclidean distance from each point in ``A`` to ``p``.
+ """
+ C = []
+ for a in A:
+ C.append(((a[0] - p[0])**2 + (a[1] - p[1])**2 + (a[2] - p[2])**2)**(1/2))
+ return C
+
+def define_patch(points,central_point,patch_size):
+ """
+ Select the nearest ``patch_size`` points around a central point.
+
+ Args:
+ points ('iterable'):
+ Iterable of candidate 3D points.
+
+ central_point ('iterable'):
+ 3D point used as the patch center.
+
+ patch_size ('int'):
+ Number of points to include in the patch.
+
+ Returns:
+ ('tuple[list[float], list[tuple[float, float, float]]]'):
+ Pair with:
+ - Distances from all input points to ``central_point``.
+ - Coordinates of the selected patch points.
+ """
+ site_positions = []
+ distance_to_central_point = calculate_distance_vector_point(points,central_point)
+ points_index = sorted(range(len(distance_to_central_point)), key=lambda sub: distance_to_central_point[sub])[:patch_size]
+ for index in points_index:
+ site_positions.append((points[index][0],points[index][1],points[index][2]))
+ return distance_to_central_point, site_positions
+
+def check_patch_overlaps(sites_positions,number_patches):
+ """
+ Check for overlapping site coordinates between patches.
+
+ Args:
+ sites_positions ('list'):
+ List containing one list of coordinates per patch.
+
+ number_patches ('int'):
+ Number of patches to compare.
+
+ Returns:
+ ('int'):
+ Returns ``0`` when no overlap is detected.
+
+ Raises:
+ ValueError:
+ If overlapping coordinates are found between any pair of patches.
+ """
+ overlapped_sites = []
+ for i in range(number_patches-1):
+ for j in range(i + 1, number_patches):
+ overlapped_sites.append(set(sites_positions[i]) & set(sites_positions[j]))
+ for overlap in overlapped_sites:
+ if len(overlap) != 0:
+ raise ValueError("The patchies are overlapping in {} sites. Please adjust the angle between them.\n".format(len(overlap)));
+ else:
+ 0
+ return 0
+
+def calculate_distance_between_points_on_sphere(points):
+ """
+ Compute nearest-neighbor distance statistics for points on a sphere.
+
+ Args:
+ points ('iterable'):
+ Nested iterable with 3D point coordinates.
+
+ Returns:
+ ('tuple[float, float, float]'):
+ Tuple ``(mean, std, stderr)`` of nearest-neighbor distances.
+ """
+ points = np.vstack(points)
+ tree = cKDTree(points)
+ distances, _ = tree.query(points, k=7) # Assuming 6 neighbors plus the point itself k = 7
+ nearest_neighbors_dist = distances[:, 1] # Removing self (distance = 0)
+ avg_dis = np.mean(nearest_neighbors_dist)
+ dev_dis = np.std(nearest_neighbors_dist)
+ err_dis = dev_dis / np.sqrt(len(nearest_neighbors_dist))
+ return avg_dis, dev_dis, err_dis
+
+def calculate_dipole_moment(charges, positions):
+ """
+ Compute dipole moment for a set of point charges.
+
+ Args:
+ charges ('iterable'):
+ Charge values.
+
+ positions ('iterable'):
+ 3D coordinates matching ``charges``.
+
+ Returns:
+ ('tuple[numpy.ndarray, float]'):
+ Dipole vector and dipole magnitude.
+ """
+ dipole_moment = np.sum(np.array(charges)[:, None] * np.array(positions), axis=0)
+ dipole_magnitude = np.linalg.norm(dipole_moment)
+ return dipole_moment, dipole_magnitude
+
+def get_nanoparticle_properties(pmb, name):
+ """
+ Return the computed properties of a nanoparticle template.
+
+ Args:
+ pmb ('pyMBE.pymbe_library'):
+ Active pyMBE object with a populated template database.
+
+ name ('str'):
+ Name of the nanoparticle template.
+
+ Returns:
+ ('dict'): Dictionary with geometric, compositional, and electrostatic
+ properties as returned by
+ :meth:`~pyMBE.storage.templates.nanoparticle.NanoparticleTemplate.calculate_nanoparticle_properties`.
+
+ Raises:
+ ValueError:
+ If no nanoparticle template with the given name exists in the database.
+ """
+ tpl = pmb.db.get_template(name=name, pmb_type="nanoparticle")
+ return tpl.calculate_nanoparticle_properties(pmb)
+
+
+def print_nanoparticle_properties(properties, name=""):
+ """
+ Print nanoparticle properties in a human-readable format.
+
+ Args:
+ properties ('dict'):
+ Dictionary of nanoparticle properties as returned by
+ :func:`get_nanoparticle_properties`.
+
+ name ('str', optional):
+ Nanoparticle name shown in the header. Defaults to empty string.
+ """
+ header = f"Properties of nanoparticle '{name}'" if name else "Nanoparticle properties"
+ print(header)
+ for key, value in properties.items():
+ print(f" {key}: {value}")
+
+
+def calculate_quadrupole_moment(charges, positions):
+ """
+ Compute quadrupole moment tensor for a set of point charges.
+
+ Args:
+ charges ('iterable'):
+ Charge values.
+
+ positions ('iterable'):
+ 3D coordinates matching ``charges``.
+
+ Returns:
+ ('tuple[numpy.ndarray, float, numpy.ndarray]'):
+ Quadrupole tensor (3x3), Frobenius norm, and eigenvalues.
+ """
+ Q = np.zeros((3, 3))
+ positions = np.array(positions)
+ for q, r in zip(charges, positions):
+ r_outer = np.outer(r, r)
+ Q += q * (3 * r_outer - np.eye(3) * np.dot(r, r))
+ quadrupole_magnitude = np.linalg.norm(Q)
+ eigenvalues, _ = np.linalg.eigh(Q)
+ return Q, quadrupole_magnitude, eigenvalues
diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py
index f2bd213b..fb522af9 100644
--- a/pyMBE/pyMBE.py
+++ b/pyMBE/pyMBE.py
@@ -36,6 +36,7 @@
from pyMBE.storage.templates.peptide import PeptideTemplate
from pyMBE.storage.templates.protein import ProteinTemplate
from pyMBE.storage.templates.hydrogel import HydrogelTemplate, HydrogelNode, HydrogelChain
+from pyMBE.storage.templates.nanoparticle import NanoparticleTemplate
from pyMBE.storage.templates.bond import BondTemplate
from pyMBE.storage.templates.lj import LJInteractionTemplate
## Instances
@@ -46,10 +47,12 @@
from pyMBE.storage.instances.protein import ProteinInstance
from pyMBE.storage.instances.bond import BondInstance
from pyMBE.storage.instances.hydrogel import HydrogelInstance
+from pyMBE.storage.instances.nanoparticle import NanoparticleInstance
## Reactions
from pyMBE.storage.reactions.reaction import Reaction, ReactionParticipant
# Utilities
import pyMBE.lib.handy_functions as hf
+import pyMBE.lib.nanoparticle_tools as np_aux
import pyMBE.storage.io as io
class pymbe_library():
@@ -388,6 +391,8 @@ def _get_label_id_map(self, pmb_type):
label="assembly_map"
elif pmb_type in self.db._molecule_like_types:
label="molecule_map"
+ elif pmb_type == "particle":
+ label="all"
else:
label=f"{pmb_type}_map"
return label
@@ -432,9 +437,9 @@ def _get_template_type(self, name, allowed_types):
registered_pmb_types_with_name = self.db._find_template_types(name=name)
filtered_types = allowed_types.intersection(registered_pmb_types_with_name)
if len(filtered_types) > 1:
- raise ValueError(f"Ambiguous template name '{name}': found {len(filtered_types)} templates in the pyMBE database. Molecule creation aborted.")
+ raise ValueError(f"Ambiguous template name '{name}': found both 'molecule' and 'peptide' templates in the pyMBE database. Molecule creation aborted.")
if len(filtered_types) == 0:
- raise ValueError(f"No {allowed_types} template found with name '{name}'. Found templates of types: {filtered_types}.")
+ raise ValueError(f"No 'molecule' or 'peptide' template found with name '{name}'. Found templates of types: {filtered_types}.")
return next(iter(filtered_types))
def _delete_particles_from_espresso(self, particle_ids, espresso_system):
@@ -488,7 +493,9 @@ def calculate_center_of_mass(self, instance_id, pmb_type, espresso_system):
axis_list = [0,1,2]
inst = self.db.get_instance(pmb_type=pmb_type,
instance_id=instance_id)
- particle_id_list = self.get_particle_id_map(object_name=inst.name)["all"]
+ id_map = self.get_particle_id_map(object_name=inst.name)
+ label = self._get_label_id_map(pmb_type=pmb_type)
+ particle_id_list = id_map[label][instance_id]
for pid in particle_id_list:
for axis in axis_list:
center_of_mass [axis] += espresso_system.part.by_id(pid).pos[axis]
@@ -673,7 +680,7 @@ def calculate_net_charge(self,espresso_system,object_name,pmb_type,dimensionless
object_name (str):
Name of the object (e.g. molecule, residue, peptide, protein).
pmb_type (str):
- Type of object to analyze. Must be molecule-like.
+ Type of object to analyze.
dimensionless (bool, optional):
If True, return charge as a pure number.
If False, return a quantity with reduced_charge units.
@@ -684,9 +691,12 @@ def calculate_net_charge(self,espresso_system,object_name,pmb_type,dimensionless
"""
id_map = self.get_particle_id_map(object_name=object_name)
label = self._get_label_id_map(pmb_type=pmb_type)
- instance_map = id_map[label]
+ if pmb_type == "particle":
+ iterable = ((pid, [pid]) for pid in id_map[label])
+ else:
+ iterable = id_map[label].items()
charges = {}
- for instance_id, particle_ids in instance_map.items():
+ for instance_id, particle_ids in iterable:
if dimensionless:
net_charge = 0.0
else:
@@ -989,8 +999,6 @@ def create_molecule(self, name, number_of_molecules, espresso_system, list_of_fi
Notes:
- This function can be used to create both molecules and peptides.
"""
- pmb_type = self._get_template_type(name=name,
- allowed_types={"molecule", "peptide"})
if number_of_molecules <= 0:
return {}
if list_of_first_residue_positions is not None:
@@ -1002,6 +1010,8 @@ def create_molecule(self, name, number_of_molecules, espresso_system, list_of_fi
if len(list_of_first_residue_positions) != number_of_molecules:
raise ValueError(f"Number of positions provided in {list_of_first_residue_positions} does not match number of molecules desired, {number_of_molecules}")
+ pmb_type = self._get_template_type(name=name,
+ allowed_types={"molecule", "peptide"})
# Generate an arbitrary random unit vector
if backbone_vector is None:
backbone_vector = self.generate_random_points_in_a_sphere(center=[0,0,0],
@@ -1099,6 +1109,189 @@ def create_molecule(self, name, number_of_molecules, espresso_system, list_of_fi
pos_index+=1
molecule_ids.append(molecule_id)
return molecule_ids
+
+ def _create_nanoparticle_sites_positions(self, nanoparticle_tpl, tolerance=1e-6, angle_between_patches=180):
+ """
+ Build per-patch site coordinates for a nanoparticle template.
+
+ Args:
+ nanoparticle_tpl ('NanoparticleTemplate'):
+ Nanoparticle template from the pyMBE database.
+
+ tolerance ('float', optional):
+ Convergence tolerance used in the point-distribution algorithm.
+
+ angle_between_patches ('float', optional):
+ Angle (in degrees) between the two primary patches when
+ ``number_of_patches_of_primary_sites == 2``.
+
+ Returns:
+ ('list[dict]'):
+ List of dictionaries with ``particle_name``, ``positions``, and
+ ``number_of_sites`` for each site patch.
+ """
+ properties = nanoparticle_tpl.calculate_nanoparticle_properties(self)
+ if properties["total_number_of_sites"] <= 0:
+ return []
+
+ core_particle_tpl = self.db.get_template(name=nanoparticle_tpl.core_particle_name,
+ pmb_type="particle")
+ core_state = self.db.get_template(name=core_particle_tpl.initial_state,
+ pmb_type="particle_state")
+ core_radius = self.get_radius_map(dimensionless=False)[core_state.es_type]
+ sites_radius = (core_radius - (0.5 * self.units("reduced_length"))).to("reduced_length").magnitude
+
+ total_number_of_sites = properties["total_number_of_sites"]
+ number_primary_patches = nanoparticle_tpl.number_of_patches_of_primary_sites
+ number_primary_sites_per_patch = properties["number_of_primary_sites_per_patch"]
+ number_secondary_sites = properties["number_of_secondary_sites"]
+
+ root_edges = np_aux.uniform_distribution_sites_on_sphere(number_of_edges=total_number_of_sites,
+ tolerance=tolerance)
+ nanoparticle_edges = np.multiply(root_edges, sites_radius)
+
+ primary_patch_positions = []
+ if number_primary_patches <= 2:
+ initial_edge = [nanoparticle_edges[0]]
+ _, sites_positions_patch = np_aux.define_patch(points=nanoparticle_edges,
+ central_point=initial_edge[0],
+ patch_size=number_primary_sites_per_patch)
+ primary_patch_positions.append(sites_positions_patch)
+
+ if number_primary_patches == 2:
+ distance_omega = (2 * sites_radius**2 - 2 * sites_radius**2 * np.cos(np.radians(angle_between_patches)))**(1 / 2)
+ distance_to_patch_1 = list(
+ np.abs(np.array(np_aux.calculate_distance_vector_point(nanoparticle_edges, initial_edge[0])) - distance_omega)
+ )
+ second_edge = nanoparticle_edges[distance_to_patch_1.index(min(distance_to_patch_1))]
+ _, sites_positions_patch = np_aux.define_patch(points=nanoparticle_edges,
+ central_point=second_edge,
+ patch_size=number_primary_sites_per_patch)
+ primary_patch_positions.append(sites_positions_patch)
+ np_aux.check_patch_overlaps(sites_positions=primary_patch_positions,
+ number_patches=number_primary_patches)
+ elif number_primary_patches > 2:
+ initial_patch_edges = np_aux.uniform_distribution_sites_on_sphere(number_of_edges=number_primary_patches,
+ tolerance=tolerance)
+ initial_patch_scaled_edges = np.multiply(initial_patch_edges, sites_radius)
+ initial_edges = []
+ for scaled_edge in initial_patch_scaled_edges:
+ comparison = np_aux.calculate_distance_vector_point(nanoparticle_edges, scaled_edge)
+ initial_edges.append(nanoparticle_edges[comparison.index(min(comparison))])
+ for patch_index in range(number_primary_patches):
+ _, sites_positions_patch = np_aux.define_patch(points=nanoparticle_edges,
+ central_point=initial_edges[patch_index],
+ patch_size=number_primary_sites_per_patch)
+ primary_patch_positions.append(sites_positions_patch)
+ np_aux.check_patch_overlaps(sites_positions=primary_patch_positions,
+ number_patches=number_primary_patches)
+
+ remaining_positions = set(map(tuple, nanoparticle_edges))
+ for patch_positions in primary_patch_positions:
+ remaining_positions = remaining_positions.difference(set(map(tuple, patch_positions)))
+
+ sites_to_create = []
+ for patch_positions in primary_patch_positions:
+ sites_to_create.append({"particle_name": nanoparticle_tpl.primary_site_particle_name,
+ "positions": [list(position) for position in patch_positions],
+ "number_of_sites": len(patch_positions)})
+ if nanoparticle_tpl.secondary_site_particle_name is not None and number_secondary_sites > 0:
+ secondary_positions = [list(position) for position in sorted(remaining_positions)]
+ sites_to_create.append({"particle_name": nanoparticle_tpl.secondary_site_particle_name,
+ "positions": secondary_positions,
+ "number_of_sites": len(secondary_positions)})
+ return sites_to_create
+
+ def create_nanoparticle(self, name, number_of_nanoparticles, espresso_system, list_core_particle_positions=None, fix=False):
+ """
+ Creates one or more nanoparticles in an ESPResSo system using a nanoparticle
+ template from the pyMBE database.
+
+ Args:
+ name ('str'):
+ Label of a nanoparticle template in the pyMBE database.
+
+ number_of_nanoparticles ('int'):
+ Number of nanoparticle instances to create.
+
+ espresso_system ('espressomd.system.System'):
+ ESPResSo system where particles are created.
+
+ list_core_particle_positions ('list', optional):
+ Nested list with one ``[x, y, z]`` position per nanoparticle core.
+ If omitted, random core positions are used.
+
+ fix ('bool', optional):
+ If ``True``, all particles of each nanoparticle are created as fixed.
+
+ Returns:
+ ('list' of 'int'):
+ List of IDs of the created nanoparticle instances.
+
+ """
+ if number_of_nanoparticles <= 0:
+ return []
+ if not self.db._has_template(name=name, pmb_type="nanoparticle"):
+ raise ValueError(f"Nanoparticle template with name '{name}' is not defined in the pyMBE database.")
+ if list_core_particle_positions is not None:
+ if len(list_core_particle_positions) != number_of_nanoparticles:
+ raise ValueError(
+ f"Number of positions ({len(list_core_particle_positions)}) does not match "
+ f"number_of_nanoparticles ({number_of_nanoparticles})."
+ )
+ for item in list_core_particle_positions:
+ if not isinstance(item, list) or len(item) != 3:
+ raise ValueError(
+ "Each core position must be a list with three coordinates [x, y, z]."
+ )
+ nanoparticle_ids = []
+ nanoparticle_tpl = self.db.get_template(name=name, pmb_type="nanoparticle")
+ site_patch_specs = self._create_nanoparticle_sites_positions(nanoparticle_tpl=nanoparticle_tpl)
+ for nanoparticle_index in range(number_of_nanoparticles):
+ nanoparticle_id = self.db._propose_instance_id(pmb_type="nanoparticle")
+ nanoparticle_ids.append(nanoparticle_id)
+ if list_core_particle_positions is None:
+ core_particle_id = self.create_particle(name=nanoparticle_tpl.core_particle_name,
+ espresso_system=espresso_system,
+ number_of_particles=1,
+ fix=fix)[0]
+ else:
+ core_particle_id = self.create_particle(name=nanoparticle_tpl.core_particle_name,
+ espresso_system=espresso_system,
+ position=[list_core_particle_positions[nanoparticle_index]],
+ number_of_particles=1,
+ fix=fix)[0]
+ self.db._update_instance(instance_id=core_particle_id,
+ pmb_type="particle",
+ attribute="molecule_id",
+ value=nanoparticle_id)
+ core_position = np.array(espresso_system.part.by_id(core_particle_id).pos)
+ patch_ids = []
+ all_sites_ids = []
+ for patch_spec in site_patch_specs:
+ if patch_spec["number_of_sites"] <= 0:
+ patch_ids.append([])
+ continue
+ translated_positions = (np.array(patch_spec["positions"]) + core_position).tolist()
+ created_ids = self.create_particle(name=patch_spec["particle_name"],
+ espresso_system=espresso_system,
+ position=translated_positions,
+ number_of_particles=patch_spec["number_of_sites"],
+ fix=fix)
+ for particle_id in created_ids:
+ self.db._update_instance(instance_id=particle_id,
+ pmb_type="particle",
+ attribute="molecule_id",
+ value=nanoparticle_id)
+ patch_ids.append(created_ids)
+ all_sites_ids.extend(created_ids)
+ self.db._register_instance(NanoparticleInstance(name=name,
+ molecule_id=nanoparticle_id))
+ if not fix:
+ self.enable_motion_of_rigid_object(instance_id=nanoparticle_id,
+ pmb_type="nanoparticle",
+ espresso_system=espresso_system)
+ return nanoparticle_ids
def create_particle(self, name, espresso_system, number_of_particles, position=None, fix=False):
"""
@@ -1197,6 +1390,7 @@ def create_protein(self, name, number_of_proteins, espresso_system, topology_dic
return
if not self.db._has_template(name=name, pmb_type="protein"):
raise ValueError(f"Protein template with name '{name}' is not defined in the pyMBE database.")
+
protein_tpl = self.db.get_template(pmb_type="protein", name=name)
box_half = espresso_system.box_l[0] / 2.0
# Create protein
@@ -1494,6 +1688,45 @@ def define_hydrogel(self, name, node_map, chain_map):
chain_map=chains)
self.db._register_template(tpl)
+ def define_nanoparticle(self, name, core_particle_name, total_number_of_sites, primary_site_particle_name, fraction_primary_sites, number_of_patches_of_primary_sites, secondary_site_particle_name=None):
+ """
+ Defines a nanoparticle template in the pyMBE database.
+
+ Args:
+ name ('str'):
+ Unique label that identifies the nanoparticle template.
+
+ core_particle_name ('str'):
+ Name of the particle template used as the nanoparticle core.
+
+ total_number_of_sites ('int'):
+ Total number of grafting/interaction sites on the nanoparticle surface.
+ The surface density is computed from this value and the core radius.
+
+ primary_site_particle_name ('str'):
+ Particle template used for the primary site type.
+
+ fraction_primary_sites ('float'):
+ Fraction of sites assigned to the primary site type.
+
+ number_of_patches_of_primary_sites ('int'):
+ Number of primary-site patches on the nanoparticle surface.
+
+ secondary_site_particle_name ('str', optional):
+ Optional particle template used for a secondary site type.
+ Defaults to None.
+
+ """
+ tpl = NanoparticleTemplate(name=name,
+ core_particle_name=core_particle_name,
+ total_number_of_sites=total_number_of_sites,
+ primary_site_particle_name=primary_site_particle_name,
+ fraction_primary_sites=fraction_primary_sites,
+ number_of_patches_of_primary_sites=number_of_patches_of_primary_sites,
+ secondary_site_particle_name=secondary_site_particle_name)
+ self.db._register_template(tpl)
+
+
def define_molecule(self, name, residue_list):
"""
Defines a molecule template in the pyMBE database.
@@ -1899,9 +2132,33 @@ def enable_motion_of_rigid_object(self, instance_id, pmb_type, espresso_system):
center_of_mass = self.calculate_center_of_mass (instance_id=instance_id,
espresso_system=espresso_system,
pmb_type=pmb_type)
+ rigid_center_name = f"{inst.name}_rigid_center"
+ if rigid_center_name not in self.db._templates["particle"]:
+ part_tpl = ParticleTemplate(name=f"{inst.name}_rigid_center",
+ sigma=PintQuantity.from_quantity(q=self.units.Quantity(0, "nm"),
+ ureg=self.units,
+ expected_dimension="length"),
+ offset=PintQuantity.from_quantity(q=self.units.Quantity(0, "nm"),
+ ureg=self.units,
+ expected_dimension="length"),
+ cutoff=PintQuantity.from_quantity(q=self.units.Quantity(0, "nm"),
+ ureg=self.units,
+ expected_dimension="length"),
+ epsilon=PintQuantity.from_quantity(q=self.units.Quantity(0, "kJ"),
+ ureg=self.units,
+ expected_dimension="energy"),
+ initial_state="None")
+ self.db._register_template(part_tpl)
rigid_object_center = espresso_system.part.add(pos=center_of_mass,
rotation=[True,True,True],
type=self.propose_unused_type())
+ part_inst = ParticleInstance(particle_id=rigid_object_center.id,
+ name=f"{inst.name}_rigid_center",
+ residue_id=instance_id if pmb_type == "residue" else None,
+ initial_state="None",
+ molecule_id=instance_id if pmb_type in self.db._molecule_like_types else None,
+ assembly_id=instance_id if pmb_type in self.db._assembly_like_types else None,)
+ self.db._register_instance(part_inst)
rigid_object_center.mass = len(particle_ids_list)
momI = 0
for pid in particle_ids_list:
diff --git a/pyMBE/storage/instances/nanoparticle.py b/pyMBE/storage/instances/nanoparticle.py
new file mode 100644
index 00000000..24a69195
--- /dev/null
+++ b/pyMBE/storage/instances/nanoparticle.py
@@ -0,0 +1,53 @@
+#
+# Copyright (C) 2026 pyMBE-dev team
+#
+# This file is part of pyMBE.
+#
+# pyMBE is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# pyMBE is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+
+from pyMBE.storage.base_type import PMBBaseModel
+from pydantic import validator
+from typing import Optional
+
+class NanoparticleInstance(PMBBaseModel):
+ """
+ Persistent instance representation of a nanoparticle.
+
+ Attributes:
+ pmb_type ('str'):
+ Fixed string identifying this object as a nanoparticle instance.
+ Always ``"nanoparticle"``.
+
+ name ('str'):
+ Name of the nanoparticle **template** from which this instance was
+ created. This must correspond to an existing
+ ``NanoparticleTemplate`` in the database.
+
+ molecule_id ('int'):
+ Unique non-negative integer identifying the nanoparticle instance within the database.
+
+ assembly_id (int | None):
+ Identifier of the super-parent assembly (e.g. hydrogel) to which this nanoparticle belongs. ``None`` indicates that the nanoparticle is not assigned to any assembly.
+ """
+ pmb_type: str = "nanoparticle"
+ name: str
+ molecule_id: int
+ assembly_id: Optional[int] = None
+
+ @validator("molecule_id")
+ def validate_molecule_id(cls, mid):
+ if mid < 0:
+ raise ValueError("molecule_id must be a non-negative integer.")
+ return mid
diff --git a/pyMBE/storage/io.py b/pyMBE/storage/io.py
index 70826528..bae94659 100644
--- a/pyMBE/storage/io.py
+++ b/pyMBE/storage/io.py
@@ -40,6 +40,8 @@
from pyMBE.storage.instances.protein import ProteinInstance
from pyMBE.storage.templates.hydrogel import HydrogelTemplate, HydrogelNode, HydrogelChain
from pyMBE.storage.instances.hydrogel import HydrogelInstance
+from pyMBE.storage.templates.nanoparticle import NanoparticleTemplate
+from pyMBE.storage.instances.nanoparticle import NanoparticleInstance
from pyMBE.storage.templates.lj import LJInteractionTemplate
def _decode(s):
@@ -123,7 +125,8 @@ def _load_database_csv(db, folder):
Notes:
- PintQuantity objects are reconstructed from their dictionary representation.
- - Supports particle, residue, molecule, peptide, protein, bond, and hydrogel types.
+ - Supports particle, residue, molecule, peptide, protein, bond, hydrogel,
+ nanoparticle, and lj types.
"""
folder = Path(folder)
if not folder.exists():
@@ -137,6 +140,7 @@ def _load_database_csv(db, folder):
"peptide",
"protein",
"hydrogel",
+ "nanoparticle",
"lj"]
# TEMPLATES
for pmb_type in pyMBE_types:
@@ -228,6 +232,16 @@ def _load_database_csv(db, folder):
node_map=node_map,
chain_map=chain_map)
templates[tpl.name] = tpl
+ elif pmb_type == "nanoparticle":
+ secondary_site = row.get("secondary_site_particle_name", "") or None
+ tpl = NanoparticleTemplate(name=row["name"],
+ core_particle_name=row["core_particle_name"],
+ total_number_of_sites=int(row["total_number_of_sites"]),
+ primary_site_particle_name=row["primary_site_particle_name"],
+ fraction_primary_sites=float(row["fraction_primary_sites"]),
+ number_of_patches_of_primary_sites=int(row["number_of_patches_of_primary_sites"]),
+ secondary_site_particle_name=secondary_site)
+ templates[tpl.name] = tpl
elif pmb_type == "lj":
sigma_d = _decode(row["sigma"])
epsilon_d = _decode(row["epsilon"])
@@ -312,6 +326,13 @@ def _load_database_csv(db, folder):
inst = HydrogelInstance(name=row["name"],
assembly_id=int(row["assembly_id"]))
instances[inst.assembly_id] = inst
+ elif pmb_type == "nanoparticle":
+ molecule_val = row.get("molecule_id", "") or ""
+ assembly_val = row.get("assembly_id", "") or ""
+ inst = NanoparticleInstance(name=row["name"],
+ molecule_id=int(molecule_val),
+ assembly_id=None if assembly_val == "" else int(assembly_val))
+ instances[inst.molecule_id] = inst
db._instances[pmb_type] = instances
# REACTIONS
@@ -406,6 +427,14 @@ def _save_database_csv(db, folder):
rows.append({"name": tpl.name,
"node_map": _encode([node.dict() for node in tpl.node_map]),
"chain_map": _encode([chain.dict() for chain in tpl.chain_map])})
+ elif pmb_type == "nanoparticle" and isinstance(tpl, NanoparticleTemplate):
+ rows.append({"name": tpl.name,
+ "core_particle_name": tpl.core_particle_name,
+ "total_number_of_sites": tpl.total_number_of_sites,
+ "primary_site_particle_name": tpl.primary_site_particle_name,
+ "fraction_primary_sites": tpl.fraction_primary_sites,
+ "number_of_patches_of_primary_sites": tpl.number_of_patches_of_primary_sites,
+ "secondary_site_particle_name": tpl.secondary_site_particle_name if tpl.secondary_site_particle_name is not None else ""})
# LJ TEMPLATE
elif pmb_type == "lj" and isinstance(tpl, LJInteractionTemplate):
rows.append({"name": tpl.name,
@@ -469,6 +498,11 @@ def _save_database_csv(db, folder):
rows.append({"pmb_type": pmb_type,
"name": inst.name,
"assembly_id": int(inst.assembly_id)})
+ elif pmb_type == "nanoparticle" and isinstance(inst, NanoparticleInstance):
+ rows.append({"pmb_type": pmb_type,
+ "name": inst.name,
+ "molecule_id": int(inst.molecule_id),
+ "assembly_id": int(inst.assembly_id) if inst.assembly_id is not None else ""})
else:
# fallback to dict
try:
@@ -488,4 +522,4 @@ def _save_database_csv(db, folder):
"reaction_type": rx.reaction_type,
"metadata": _encode(rx.metadata) if getattr(rx, "metadata", None) is not None else ""})
if rows:
- pd.DataFrame(rows).to_csv(os.path.join(folder, "reactions.csv"), index=False)
\ No newline at end of file
+ pd.DataFrame(rows).to_csv(os.path.join(folder, "reactions.csv"), index=False)
diff --git a/pyMBE/storage/manager.py b/pyMBE/storage/manager.py
index f4b3f404..ecbffc4a 100644
--- a/pyMBE/storage/manager.py
+++ b/pyMBE/storage/manager.py
@@ -35,6 +35,8 @@
from pyMBE.storage.instances.protein import ProteinInstance
from pyMBE.storage.templates.hydrogel import HydrogelTemplate
from pyMBE.storage.instances.hydrogel import HydrogelInstance
+from pyMBE.storage.templates.nanoparticle import NanoparticleTemplate
+from pyMBE.storage.instances.nanoparticle import NanoparticleInstance
from pyMBE.storage.templates.lj import LJInteractionTemplate
from pyMBE.storage.pint_quantity import PintQuantity
@@ -86,7 +88,8 @@ def __init__(self,units):
self._reactions: Dict[str, Reaction] = {}
self._molecule_like_types = ["molecule",
"peptide",
- "protein"]
+ "protein",
+ "nanoparticle"]
self._assembly_like_types = ["hydrogel"]
self._pmb_types = ["particle", "residue"] + self._molecule_like_types + self._assembly_like_types
self.espresso_bond_instances= {}
@@ -146,6 +149,12 @@ def _collect_particle_templates(self, name, pmb_type):
if pmb_type in self._molecule_like_types:
tpl = self.get_template(name=name,
pmb_type=pmb_type)
+ if pmb_type == "nanoparticle":
+ counts[tpl.core_particle_name] += 1
+ counts[tpl.primary_site_particle_name] += 1
+ if tpl.secondary_site_particle_name is not None:
+ counts[tpl.secondary_site_particle_name] += 1
+ return counts
for res_name in tpl.residue_list:
sub = self._collect_particle_templates(name=res_name,
pmb_type="residue")
@@ -355,6 +364,15 @@ def _get_templates_df(self, pmb_type):
"particle_name1": tpl.particle_name1,
"particle_name2": tpl.particle_name2,
"parameters": parameters})
+ elif pmb_type == "nanoparticle":
+ rows.append({"pmb_type": tpl.pmb_type,
+ "name": tpl.name,
+ "core_particle_name": tpl.core_particle_name,
+ "total_number_of_sites": tpl.total_number_of_sites,
+ "primary_site_particle_name": tpl.primary_site_particle_name,
+ "fraction_primary_sites": tpl.fraction_primary_sites,
+ "number_of_patches_of_primary_sites": tpl.number_of_patches_of_primary_sites,
+ "secondary_site_particle_name": tpl.secondary_site_particle_name})
else:
# Generic representation for other types
rows.append(tpl.dict())
@@ -431,6 +449,9 @@ def _register_instance(self, instance):
elif isinstance(instance, HydrogelInstance):
pmb_type = "hydrogel"
iid = instance.assembly_id
+ elif isinstance(instance, NanoparticleInstance):
+ pmb_type = "nanoparticle"
+ iid = instance.molecule_id
else:
raise TypeError("Unsupported instance type")
self._instances.setdefault(pmb_type, {})
@@ -612,7 +633,7 @@ def _propose_instance_id(self, pmb_type):
if not used_ids:
return 0
else:
- if pmb_type not in self._instances:
+ if pmb_type not in self._instances or len(self._instances[pmb_type]) == 0:
return 0
used_ids = list(self._instances[pmb_type].keys())
return max(used_ids) + 1
diff --git a/pyMBE/storage/pint_quantity.py b/pyMBE/storage/pint_quantity.py
index 3491dc8b..e3e576d1 100644
--- a/pyMBE/storage/pint_quantity.py
+++ b/pyMBE/storage/pint_quantity.py
@@ -21,9 +21,10 @@
import pint
# dimension -> representative unit used to check dimensionality
_DIMENSION_REPRESENTATIVE = {"length": "nm",
- "energy": "meV",
- "energy/length**2": "meV/nm**2",
- "dimensionless": "dimensionless",} # extend as needed
+ "energy": "meV",
+ "energy/length**2": "meV/nm**2",
+ "length**-2": "nm**-2",
+ "dimensionless": "dimensionless",} # extend as needed
@dataclass
class PintQuantity:
diff --git a/pyMBE/storage/templates/nanoparticle.py b/pyMBE/storage/templates/nanoparticle.py
new file mode 100644
index 00000000..ae03ac48
--- /dev/null
+++ b/pyMBE/storage/templates/nanoparticle.py
@@ -0,0 +1,181 @@
+#
+# Copyright (C) 2026 pyMBE-dev team
+#
+# This file is part of pyMBE.
+#
+# pyMBE is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# pyMBE is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+
+from pyMBE.storage.base_type import PMBBaseModel
+from pydantic import Field
+import numpy as np
+
+class NanoparticleTemplate(PMBBaseModel):
+ """
+ Template defining a nanoparticle in the pyMBE database.
+
+ Attributes:
+ pmb_type ('str'):
+ Fixed type identifier. Always "nanoparticle".
+
+ name ('str'):
+ Unique name of the nanoparticle template.
+
+ core_particle_name ('str'):
+ Name of the particle template used as the nanoparticle core.
+
+ total_number_of_sites ('int'):
+ Total number of grafting/interaction sites on the nanoparticle surface.
+ The surface density is computed from this value and the core radius.
+
+ primary_site_particle_name ('str'):
+ Name of the particle template used for the primary site type.
+
+ fraction_primary_sites ('float'):
+ Fraction of all surface sites assigned to the primary site type.
+ Expected range is typically between 0 and 1.
+
+ number_of_patches_of_primary_sites ('int'):
+ Number of surface patches that contain the primary site type.
+
+ secondary_site_particle_name ('str | None'):
+ Optional particle template name for a secondary site type.
+ If not provided, only a single site type is used.
+
+ """
+ pmb_type: str = Field(default="nanoparticle", frozen=True)
+ name: str
+ core_particle_name: str
+ total_number_of_sites: int
+ primary_site_particle_name: str
+ fraction_primary_sites: float
+ number_of_patches_of_primary_sites: int
+ secondary_site_particle_name: str | None = None
+
+ def calculate_nanoparticle_properties(self, pmb):
+ """
+ Compute various nanoparticle properties from template parameters.
+
+ The method uses:
+ - Core radius from the core particle template.
+ - Site counts inferred from surface area and site surface density.
+ - Site charges taken from each site's particle initial state.
+
+ Args:
+ pmb ('pyMBE.pymbe_library'):
+ Active pyMBE object with a populated template database and unit registry.
+
+ Returns:
+ ('dict'): Dictionary with geometric, compositional, and electrostatic properties:
+ - ``nanoparticle_surface_area`` ('pint.Quantity')
+ - ``nanoparticle_volume`` ('pint.Quantity')
+ - ``total_number_of_sites`` ('int')
+ - ``real_surface_density_of_sites`` ('pint.Quantity')
+ - ``number_of_primary_sites`` ('int')
+ - ``number_of_primary_sites_per_patch`` ('int')
+ - ``number_of_secondary_sites`` ('int')
+ - ``real_fraction_primary_sites`` ('float')
+ - ``primary_site_charge_number`` ('int')
+ - ``secondary_site_charge_number`` ('int')
+ - ``total_charge`` ('pint.Quantity')
+ - ``surface_charge_density`` ('pint.Quantity')
+ - ``volume_charge_density`` ('pint.Quantity')
+
+ Notes:
+ - If ``secondary_site_particle_name`` is not set, all sites are assigned
+ to the primary site type.
+ - Primary-site counts are rounded to ensure an integer number of sites
+ per patch and exact patch occupancy.
+ """
+ if not (0.0 <= self.fraction_primary_sites <= 1.0):
+ raise ValueError("fraction_primary_sites must be between 0 and 1.")
+ if self.number_of_patches_of_primary_sites <= 0:
+ raise ValueError("number_of_patches_of_primary_sites must be > 0.")
+
+ def _get_initial_state_charge_number(particle_name: str) -> int:
+ particle_tpl = pmb.db.get_template(name=particle_name, pmb_type="particle")
+ state_name = particle_tpl.initial_state
+ if state_name is None:
+ particle_states = pmb.db.get_particle_states_templates(particle_name=particle_name)
+ if not particle_states:
+ raise ValueError(f"Particle '{particle_name}' has no defined states.")
+ state_name = next(iter(particle_states.keys()))
+ state_tpl = pmb.db.get_template(name=state_name, pmb_type="particle_state")
+ return state_tpl.z
+
+ core_particle_tpl = pmb.db.get_template(name=self.core_particle_name, pmb_type="particle")
+ core_initial_state_name = core_particle_tpl.initial_state
+ if core_initial_state_name is None:
+ raise ValueError(
+ f"Core particle '{self.core_particle_name}' has no initial_state set."
+ )
+ core_initial_state = pmb.db.get_template(name=core_initial_state_name, pmb_type="particle_state")
+ core_radius = pmb.get_radius_map(dimensionless=False)[core_initial_state.es_type]
+
+ nanoparticle_surface_area = 4.0 * np.pi * core_radius**2
+ nanoparticle_volume = (4.0 / 3.0) * np.pi * core_radius**3
+
+ total_number_of_sites = self.total_number_of_sites
+ real_surface_density_of_sites = total_number_of_sites / nanoparticle_surface_area
+
+ if self.secondary_site_particle_name is None:
+ number_of_primary_sites = total_number_of_sites
+ number_of_primary_sites_per_patch = int(
+ np.round(number_of_primary_sites / self.number_of_patches_of_primary_sites)
+ )
+ real_number_of_primary_sites = (
+ number_of_primary_sites_per_patch * self.number_of_patches_of_primary_sites
+ )
+ number_of_secondary_sites = 0
+ secondary_site_charge_number = 0
+ else:
+ number_of_primary_sites = int(np.round(total_number_of_sites * self.fraction_primary_sites))
+ number_of_primary_sites_per_patch = int(
+ np.round(number_of_primary_sites / self.number_of_patches_of_primary_sites)
+ )
+ real_number_of_primary_sites = (
+ number_of_primary_sites_per_patch * self.number_of_patches_of_primary_sites
+ )
+ number_of_secondary_sites = total_number_of_sites - real_number_of_primary_sites
+ secondary_site_charge_number = _get_initial_state_charge_number(
+ self.secondary_site_particle_name
+ )
+
+ if total_number_of_sites > 0:
+ real_fraction_primary_sites = real_number_of_primary_sites / total_number_of_sites
+ else:
+ real_fraction_primary_sites = 0.0
+
+ primary_site_charge_number = _get_initial_state_charge_number(self.primary_site_particle_name)
+ total_charge_number = (
+ real_number_of_primary_sites * primary_site_charge_number
+ + number_of_secondary_sites * secondary_site_charge_number
+ )
+ total_charge = total_charge_number * pmb.units("reduced_charge")
+ surface_charge_density = total_charge / nanoparticle_surface_area
+ volume_charge_density = total_charge / nanoparticle_volume
+
+ return {"nanoparticle_surface_area": nanoparticle_surface_area,
+ "nanoparticle_volume": nanoparticle_volume,
+ "total_number_of_sites": total_number_of_sites,
+ "real_surface_density_of_sites": real_surface_density_of_sites,
+ "number_of_primary_sites": real_number_of_primary_sites,
+ "number_of_primary_sites_per_patch": number_of_primary_sites_per_patch,
+ "number_of_secondary_sites": number_of_secondary_sites,
+ "real_fraction_primary_sites": real_fraction_primary_sites,
+ "primary_site_charge_number": primary_site_charge_number,
+ "secondary_site_charge_number": secondary_site_charge_number,
+ "total_charge": total_charge,
+ "surface_charge_density": surface_charge_density,
+ "volume_charge_density": volume_charge_density,}
diff --git a/samples/nanoparticles_grxmc.py b/samples/nanoparticles_grxmc.py
new file mode 100644
index 00000000..60c59e86
--- /dev/null
+++ b/samples/nanoparticles_grxmc.py
@@ -0,0 +1,359 @@
+#
+# Copyright (C) 2026 pyMBE-dev team
+#
+# This file is part of pyMBE.
+#
+# pyMBE is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# pyMBE is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+#Load espresso, pyMBE and other necessary libraries
+import espressomd
+from pathlib import Path
+import pandas as pd
+import argparse
+from espressomd.io.writer import vtf
+import pyMBE
+from pyMBE.lib.analysis import built_output_name
+from pyMBE.lib.handy_functions import do_reaction, setup_electrostatic_interactions, relax_espresso_system, generate_lattice_positions
+from pyMBE.lib.nanoparticle_tools import get_nanoparticle_properties, print_nanoparticle_properties
+
+# Create an instance of pyMBE library
+pmb = pyMBE.pymbe_library(seed=42)
+
+# Command line arguments
+
+parser = argparse.ArgumentParser(description='Script that runs a simulation of an ideal peptide mixture in the grand-reaction ensemble using pyMBE and ESPResSo.')
+parser.add_argument('--mode',
+ type=str,
+ default= "standard",
+ choices=["standard", "unified"],
+ help='Set if the grand-reaction method is used with unified ions or not')
+parser.add_argument('--test',
+ default=False,
+ action='store_true',
+ help='to run a short simulation for testing the script')
+parser.add_argument('--pH',
+ type=float,
+ default=7,
+ help='pH of the solution')
+parser.add_argument('--output',
+ type=Path,
+ required= False,
+ default=Path(__file__).parent / "time_series" / "nanoparticle_grxmc",
+ help='output directory')
+parser.add_argument('--no_verbose', action='store_false', help="Switch to deactivate verbose",default=True)
+args = parser.parse_args()
+
+# The trajectories of the simulations will be stored using espresso built-up functions in separed files in the folder 'frames'
+frames_path = args.output / "frames"
+frames_path.mkdir(parents=True, exist_ok=True)
+
+# Simulation parameters
+verbose = args.no_verbose
+pmb.set_reduced_units(unit_length=0.4*pmb.units.nm)
+N_samples = 1000 # to make the demonstration quick, we set this to a very low value
+MD_steps_per_sample = 1000
+N_samples_print = 1 # Write the trajectory every 100 samples
+langevin_seed = 42
+dt = 0.001
+solvent_permitivity = 78.3
+pH_value= args.pH
+ideal = False # Set to True to not consider electrostatic interactions in the system, and only sample the reactions
+
+# Nanoparticle parameters
+vol_frac_of_nanoparticles = 0.1 # Volume fraction of the nanoparticle
+number_of_nanoparticles = 20 # Total number of the nanoparticles
+nanoparticle_diameter = 4*pmb.units.reduced_length # Diameter of the nanoparticle in reduced units
+total_number_of_sites = 10 # Equivalent to 0.2 sites/reduced_length^2 for a diameter-4 nanoparticle
+pka_A_site = 4.0
+pka_B_site = 10.0
+nanoparticle_lattice_type = "fcc"
+
+# Names for the componentes of the nanoparticles
+core_particle = "core_particle"
+A_site = "A_site"
+B_site = "B_site"
+
+# Patchy distribution of sites A and B
+sites_distribution = {"main" : {"particle_name" : A_site,
+ "fraction" : 0.5,
+ "number_of_patches" : 2},
+ "secondary": {"particle_name" : B_site}}
+
+# LJ parameters for the nanoparticles
+sigma_core_particle = 1*pmb.units('reduced_length')
+sigma_sites = 0*pmb.units('reduced_length')
+epsilon = 1*pmb.units('reduced_energy')
+offset_core_particle = nanoparticle_diameter-sigma_core_particle
+cutoff_core_particle = 2**(1/6)*sigma_core_particle
+
+# Short simulation setup for testing
+
+if args.test:
+ MD_steps_per_sample = 100
+ N_samples = 10
+ phi_np = 0.1
+ np_diameter = 4
+
+# Defines the components of the nanoparticle (core particle, A and B type of sites) in the pyMBE data frame
+
+pmb.define_particle(name = core_particle,
+ z = 0,
+ sigma = sigma_core_particle,
+ epsilon = epsilon,
+ offset = offset_core_particle,
+ cutoff = cutoff_core_particle)
+
+pmb.define_particle(name = A_site,
+ acidity = "acidic",
+ pka = pka_A_site,
+ sigma = sigma_sites,
+ epsilon = epsilon)
+
+pmb.define_particle(name = B_site,
+ acidity = "basic",
+ pka = pka_B_site,
+ sigma = sigma_sites,
+ epsilon = epsilon)
+
+nanoparticle_name = "nanoparticle"
+pmb.define_nanoparticle(name = nanoparticle_name,
+ core_particle_name = core_particle,
+ total_number_of_sites = total_number_of_sites,
+ primary_site_particle_name = A_site,
+ fraction_primary_sites = sites_distribution["main"]["fraction"],
+ number_of_patches_of_primary_sites = sites_distribution["main"]["number_of_patches"],
+ secondary_site_particle_name = B_site)
+
+# Saline solution parameters
+
+c_salt = 5e-3 * pmb.units.mol/ pmb.units.L
+
+if args.mode == 'standard':
+ proton_name = 'Hplus'
+ hydroxide_name = 'OHminus'
+ sodium_name = 'Na'
+ chloride_name = 'Cl'
+
+ pmb.define_particle(name = proton_name,
+ z = 1,
+ sigma = 0.35*pmb.units.nm,
+ epsilon = 1*pmb.units('reduced_energy'))
+ pmb.define_particle(name = hydroxide_name,
+ z = -1,
+ sigma = 0.35*pmb.units.nm,
+ epsilon = 1*pmb.units('reduced_energy'))
+ pmb.define_particle(name = sodium_name,
+ z = 1,
+ sigma = 0.35*pmb.units.nm,
+ epsilon = 1*pmb.units('reduced_energy'))
+ pmb.define_particle(name = chloride_name,
+ z = -1,
+ sigma = 0.35*pmb.units.nm,
+ epsilon = 1*pmb.units('reduced_energy'))
+
+elif args.mode == 'unified':
+ cation_name = 'Na'
+ anion_name = 'Cl'
+
+ pmb.define_particle(name = cation_name,
+ z = 1,
+ sigma = 0.35*pmb.units.nm,
+ epsilon = 1*pmb.units('reduced_energy'))
+ pmb.define_particle(name = anion_name,
+ z = -1,
+ sigma = 0.35*pmb.units.nm,
+ epsilon = 1*pmb.units('reduced_energy'))
+
+# System parameters
+properties = get_nanoparticle_properties(pmb, nanoparticle_name)
+nanoparticle_volume = properties["nanoparticle_volume"].to(pmb.units('reduced_length**3'))
+volume = number_of_nanoparticles * nanoparticle_volume / vol_frac_of_nanoparticles
+L = volume ** (1./3.) # Length of the simulation box
+
+# Print nanoparticle properties
+if verbose:
+ print_nanoparticle_properties(properties, name=nanoparticle_name)
+
+
+# Create an instance of an espresso system
+
+espresso_system = espressomd.System (box_l = [L.to('reduced_length').magnitude]*3)
+
+# Create non-overlapping nanoparticle core positions on a lattice
+nanoparticle_positions = generate_lattice_positions(
+ lattice_type=nanoparticle_lattice_type,
+ number_of_sites=number_of_nanoparticles,
+ box_length=L.to('reduced_length').magnitude,
+)
+
+nanoparticle_ids = pmb.create_nanoparticle(name=nanoparticle_name,
+ espresso_system=espresso_system,
+ number_of_nanoparticles=number_of_nanoparticles,
+ list_core_particle_positions=nanoparticle_positions)
+
+if args.mode == 'standard':
+ pmb.create_counterions(object_name=nanoparticle_name,
+ cation_name=proton_name,
+ anion_name=hydroxide_name,
+ espresso_system=espresso_system) # Create counterions for the peptide chains with sequence 1
+ c_salt_calculated = pmb.create_added_salt(espresso_system=espresso_system,
+ cation_name=sodium_name,
+ anion_name=chloride_name,
+ c_salt=c_salt)
+elif args.mode == 'unified':
+ pmb.create_counterions(object_name=nanoparticle_name,
+ cation_name=cation_name,
+ anion_name=anion_name,
+ espresso_system=espresso_system) # Create counterions for the peptide chains with sequence 1
+ c_salt_calculated = pmb.create_added_salt(espresso_system=espresso_system,
+ cation_name=cation_name,
+ anion_name=anion_name,
+ c_salt=c_salt)
+
+with open(frames_path / "trajectory0.vtf", mode='w+t') as coordinates:
+ vtf.writevsf(espresso_system, coordinates)
+ vtf.writevcf(espresso_system, coordinates)
+
+# count acid/base particles
+pka_set = pmb.get_pka_set()
+acid_base_ids = []
+for name in pka_set.keys():
+ acid_base_ids+=pmb.db.find_instance_ids_by_name(pmb_type="particle",
+ name=name)
+total_ionisable_groups = len(acid_base_ids)
+
+# Get nanoparticle net charge
+if verbose:
+ print("The box length of your system is", L.to('reduced_length'), L.to('nm'))
+
+if args.mode == 'standard':
+ grxmc, ionic_strength_res = pmb.setup_grxmc_reactions(pH_res=pH_value,
+ c_salt_res=c_salt,
+ proton_name=proton_name,
+ hydroxide_name=hydroxide_name,
+ salt_cation_name=sodium_name,
+ salt_anion_name=chloride_name,
+ activity_coefficient=lambda x: 1.0,
+ use_exclusion_radius_per_type=True)
+elif args.mode == 'unified':
+ grxmc, ionic_strength_res = pmb.setup_grxmc_unified(pH_res=pH_value,
+ c_salt_res=c_salt,
+ cation_name=cation_name,
+ anion_name=anion_name,
+ activity_coefficient=lambda x: 1.0,
+ use_exclusion_radius_per_type=True)
+
+if verbose:
+ print(pmb.get_reactions_df())
+
+# Setup espresso to track the ionization of the acid/basic groups in nanoparticle sites
+type_map =pmb.get_type_map()
+types = list (type_map.values())
+espresso_system.setup_type_map(type_list = types)
+
+# Setup the non-interacting type for speeding up the sampling of the reactions
+non_interacting_type = max(type_map.values())+1
+grxmc.set_non_interacting_type (type=non_interacting_type)
+if verbose:
+ print('The non interacting type is set to ', non_interacting_type)
+
+espresso_system.time_step = dt
+
+#Save the initial state
+with open(frames_path / "trajectory1.vtf", mode='w+t') as coordinates:
+ vtf.writevsf(espresso_system, coordinates)
+ vtf.writevcf(espresso_system, coordinates)
+
+# Setup espresso to do langevin dynamics
+espresso_system.time_step= dt
+espresso_system.integrator.set_vv()
+espresso_system.thermostat.set_langevin(kT=pmb.kT.to('reduced_energy').magnitude, gamma=0.1, seed=langevin_seed)
+espresso_system.cell_system.skin=0.4
+if not ideal:
+ ##Setup the potential energy
+ if verbose:
+ print('Setup LJ interaction (this can take a few seconds)')
+ pmb.setup_lj_interactions (espresso_system=espresso_system)
+ if verbose:
+ print('Minimize energy before adding electrostatics')
+ relax_espresso_system(espresso_system=espresso_system,
+ seed=langevin_seed)
+ if verbose:
+ print('Setup and tune electrostatics (this can take a few seconds)')
+ setup_electrostatic_interactions(units=pmb.units,
+ espresso_system=espresso_system,
+ kT=pmb.kT,
+ verbose=verbose)
+ if verbose:
+ print('Minimize energy after adding electrostatics')
+ relax_espresso_system(espresso_system=espresso_system,
+ seed=langevin_seed)
+
+#Save the pyMBE dataframe in a CSV file
+#Save the pyMBE database
+pmb.save_database(folder=args.output / 'database')
+
+time_series={}
+for label in ["time", "net_charge_nanoparticle", "mean_charge_primary_sites","mean_charge_secondary_sites", "num_plus","xi_plus"]:
+ time_series[label]=[]
+
+# Main simulation loop
+N_frame=0
+for step in range(N_samples):
+ print(f"Sample {step+1}/{N_samples}")
+ if not ideal:
+ espresso_system.integrator.run(steps=MD_steps_per_sample)
+ do_reaction(grxmc, steps=total_ionisable_groups)
+ time_series["time"].append(espresso_system.time)
+ # Get net charge of nanoparticle and peptide2
+ charge_dict_nanoparticle=pmb.calculate_net_charge(espresso_system=espresso_system,
+ object_name=nanoparticle_name,
+ pmb_type="nanoparticle",
+ dimensionless=True)
+ charge_dict_A_site=pmb.calculate_net_charge(espresso_system=espresso_system,
+ object_name=A_site,
+ pmb_type="particle",
+ dimensionless=True)
+ charge_dict_B_site=pmb.calculate_net_charge(espresso_system=espresso_system,
+ object_name=B_site,
+ pmb_type="particle",
+ dimensionless=True)
+ time_series["net_charge_nanoparticle"].append(charge_dict_nanoparticle["mean"])
+ time_series["mean_charge_primary_sites"].append(charge_dict_A_site["mean"])
+ time_series["mean_charge_secondary_sites"].append(charge_dict_B_site["mean"])
+ # Get degree of ionization of primary and secondary sites
+ if args.mode == 'standard':
+ num_plus = espresso_system.number_of_particles(type=type_map["Na"])+espresso_system.number_of_particles(type=type_map["Hplus"])
+ elif args.mode == 'unified':
+ num_plus = espresso_system.number_of_particles(type=type_map["Na"])
+ time_series["num_plus"].append(num_plus)
+ concentration_plus = (num_plus/(pmb.N_A * L**3)).to("mol/L")
+ xi_plus = (concentration_plus/ionic_strength_res).magnitude
+ time_series["xi_plus"].append(xi_plus)
+ if step % N_samples_print == 0:
+ N_frame+=1
+ with open(frames_path / f"trajectory{N_frame}.vtf", mode='w+t') as coordinates:
+ vtf.writevsf(espresso_system, coordinates)
+ vtf.writevcf(espresso_system, coordinates)
+
+# Store time series
+data_path=args.output
+data_path.mkdir(parents=True, exist_ok=True)
+time_series=pd.DataFrame(time_series)
+
+filename=built_output_name(input_dict={"mode":args.mode,
+ "pH":args.pH})
+
+time_series.to_csv(data_path / f"{filename}_time_series.csv",
+ index=False)
diff --git a/testsuite/CTestTestfile.cmake b/testsuite/CTestTestfile.cmake
index aa66fb39..0f958add 100644
--- a/testsuite/CTestTestfile.cmake
+++ b/testsuite/CTestTestfile.cmake
@@ -76,3 +76,4 @@ pymbe_add_test(PATH globular_protein_unit_tests.py)
pymbe_add_test(PATH lattice_builder.py)
pymbe_add_test(PATH hydrogel_builder.py)
pymbe_add_test(PATH database_unit_tests.py)
+pymbe_add_test(PATH nanoparticle_unit_tests.py)
diff --git a/testsuite/nanoparticle_unit_tests.py b/testsuite/nanoparticle_unit_tests.py
new file mode 100644
index 00000000..fd486d3d
--- /dev/null
+++ b/testsuite/nanoparticle_unit_tests.py
@@ -0,0 +1,539 @@
+#
+# Copyright (C) 2026 pyMBE-dev team
+#
+# This file is part of pyMBE.
+#
+# pyMBE is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# pyMBE is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+
+import unittest as ut
+import tempfile
+
+import espressomd
+import pyMBE
+import pyMBE.lib.nanoparticle_tools as nanoparticle_tools
+from pyMBE.storage.instances.nanoparticle import NanoparticleInstance
+from pyMBE.storage.templates.particle import ParticleTemplate
+from pyMBE.storage.pint_quantity import PintQuantity
+
+
+# ESPResSo allows only one System instance per process.
+espresso_system = espressomd.System(box_l=[40, 40, 40])
+
+
+class TestNanoparticleCreation(ut.TestCase):
+
+ def setUp(self):
+ """
+ Reset ESPResSo state before each unit test.
+
+ Notes:
+ - ESPResSo allows one ``System`` instance per process in this test
+ file, so particles are cleared between tests.
+ """
+ espresso_system.part.clear()
+
+ def _build_pmb_with_particles(self):
+ """
+ Build a pyMBE object with particle templates used by nanoparticle tests.
+
+ Returns:
+ ('pyMBE.pymbe_library'):
+ Configured pyMBE object containing templates for ``core``,
+ ``A``, and ``B`` particles.
+ """
+ pmb = pyMBE.pymbe_library(seed=42)
+ pmb.set_reduced_units(unit_length=0.4 * pmb.units.nm,
+ Kw=1e-14)
+
+ pmb.define_particle(name="core",
+ z=0,
+ sigma=4.0 * pmb.units("reduced_length"),
+ epsilon=1.0 * pmb.units("reduced_energy"))
+ pmb.define_particle(name="A",
+ z=-1,
+ sigma=1.0 * pmb.units("reduced_length"),
+ epsilon=1.0 * pmb.units("reduced_energy"))
+ pmb.define_particle(name="B",
+ z=1,
+ sigma=1.0 * pmb.units("reduced_length"),
+ epsilon=1.0 * pmb.units("reduced_energy"))
+ return pmb
+
+ def test_create_nanoparticle_registers_instances(self):
+ """
+ Unit test: verify nanoparticle creation and database bookkeeping.
+
+ Notes:
+ - Checks nanoparticle instance IDs, particle-to-molecule linkage, and
+ particle count consistency with template-derived properties.
+ """
+ pmb = self._build_pmb_with_particles()
+ nanoparticle_name = "np"
+ pmb.define_nanoparticle(name=nanoparticle_name,
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+
+ core_positions = [[5.0, 5.0, 5.0],
+ [25.0, 25.0, 25.0]]
+ created_np_ids = pmb.create_nanoparticle(name=nanoparticle_name,
+ number_of_nanoparticles=2,
+ espresso_system=espresso_system,
+ list_core_particle_positions=core_positions,
+ fix=True)
+
+ self.assertEqual(created_np_ids, [0, 1])
+
+ np_instances = pmb.db.get_instances(pmb_type="nanoparticle")
+ self.assertEqual(sorted(np_instances.keys()), [0, 1])
+ self.assertEqual(np_instances[0].molecule_id, 0)
+ self.assertEqual(np_instances[1].molecule_id, 1)
+
+ tpl = pmb.db.get_template(pmb_type="nanoparticle",
+ name=nanoparticle_name)
+ properties = tpl.calculate_nanoparticle_properties(pmb)
+ expected_particles_per_np = 1 + properties["number_of_primary_sites"] + properties["number_of_secondary_sites"]
+
+ for nanoparticle_id, expected_core_pos in enumerate(core_positions):
+ particles_in_np = pmb.db._find_instance_ids_by_attribute(pmb_type="particle",
+ attribute="molecule_id",
+ value=nanoparticle_id)
+ self.assertEqual(len(particles_in_np), expected_particles_per_np)
+
+ core_candidates = [pid for pid in particles_in_np
+ if pmb.db.get_instance(pmb_type="particle", instance_id=pid).name == "core"]
+ self.assertEqual(len(core_candidates), 1)
+ core_pid = core_candidates[0]
+
+ core_pos = list(espresso_system.part.by_id(core_pid).pos)
+ self.assertListEqual(core_pos, expected_core_pos)
+ self.assertListEqual(list(espresso_system.part.by_id(core_pid).fix), [True, True, True])
+
+ for pid in particles_in_np:
+ self.assertEqual(pmb.db.get_instance(pmb_type="particle", instance_id=pid).molecule_id,
+ nanoparticle_id)
+ self.assertListEqual(list(espresso_system.part.by_id(pid).fix), [True, True, True])
+
+ particle_id_map = pmb.get_particle_id_map(object_name=nanoparticle_name)
+ self.assertEqual(len(particle_id_map["all"]), expected_particles_per_np * 2)
+
+ def test_create_nanoparticle_input_validation_and_empty(self):
+ """
+ Unit test: verify input validation and empty-creation behavior.
+
+ Notes:
+ - Covers zero requested nanoparticles and invalid core-position inputs.
+ """
+ pmb = self._build_pmb_with_particles()
+ pmb.define_nanoparticle(name="np",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+
+ self.assertEqual(pmb.create_nanoparticle(name="np",
+ number_of_nanoparticles=0,
+ espresso_system=espresso_system),
+ [])
+
+ with self.assertRaises(ValueError):
+ pmb.create_nanoparticle(name="np",
+ number_of_nanoparticles=2,
+ espresso_system=espresso_system,
+ list_core_particle_positions=[[1.0, 2.0, 3.0]])
+
+ with self.assertRaises(ValueError):
+ pmb.create_nanoparticle(name="np",
+ number_of_nanoparticles=1,
+ espresso_system=espresso_system,
+ list_core_particle_positions=[[1.0, 2.0]])
+
+ def test_create_nanoparticle_calls_enable_motion_when_not_fixed(self):
+ """
+ Unit test: verify rigid-body motion hook is called when ``fix=False``.
+ """
+ pmb = self._build_pmb_with_particles()
+ pmb.define_nanoparticle(name="np",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+
+ calls = []
+
+ def fake_enable_motion_of_rigid_object(instance_id, pmb_type, espresso_system):
+ _ = espresso_system
+ calls.append((instance_id, pmb_type))
+
+ pmb.enable_motion_of_rigid_object = fake_enable_motion_of_rigid_object
+
+ created_np_ids = pmb.create_nanoparticle(name="np",
+ number_of_nanoparticles=1,
+ espresso_system=espresso_system,
+ fix=False)
+
+ self.assertEqual(created_np_ids, [0])
+ self.assertEqual(calls, [(0, "nanoparticle")])
+
+ def test_create_nanoparticle_sites_positions_variants(self):
+ """
+ Unit test: verify site-position generation across template variants.
+
+ Notes:
+ - Covers two-patch, multi-patch, and zero-site configurations.
+ """
+ pmb = self._build_pmb_with_particles()
+
+ # Two primary patches + secondary sites
+ pmb.define_nanoparticle(name="np_two",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+ tpl_two = pmb.db.get_template(pmb_type="nanoparticle", name="np_two")
+ properties_two = tpl_two.calculate_nanoparticle_properties(pmb)
+ sites_two = pmb._create_nanoparticle_sites_positions(nanoparticle_tpl=tpl_two)
+ self.assertEqual(len(sites_two), 3)
+ self.assertEqual(sum(item["number_of_sites"] for item in sites_two), properties_two["total_number_of_sites"])
+
+ # More than two primary patches, no secondary type.
+ # Overlap detection can be sensitive to geometric degeneracy for small systems;
+ # here we stub overlap validation to exercise this branch deterministically.
+ pmb.define_nanoparticle(name="np_three",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=1.0,
+ number_of_patches_of_primary_sites=3,
+ secondary_site_particle_name=None)
+ tpl_three = pmb.db.get_template(pmb_type="nanoparticle", name="np_three")
+ original_check_patch_overlaps = nanoparticle_tools.check_patch_overlaps
+ nanoparticle_tools.check_patch_overlaps = lambda sites_positions, number_patches: 0
+ try:
+ sites_three = pmb._create_nanoparticle_sites_positions(nanoparticle_tpl=tpl_three)
+ finally:
+ nanoparticle_tools.check_patch_overlaps = original_check_patch_overlaps
+ self.assertEqual(len(sites_three), 3)
+ self.assertTrue(all(item["particle_name"] == "A" for item in sites_three))
+
+ # Zero sites (density = 0) should return an empty specification
+ pmb.define_nanoparticle(name="np_zero",
+ core_particle_name="core",
+ total_number_of_sites=0,
+ primary_site_particle_name="A",
+ fraction_primary_sites=1.0,
+ number_of_patches_of_primary_sites=1,
+ secondary_site_particle_name=None)
+ tpl_zero = pmb.db.get_template(pmb_type="nanoparticle", name="np_zero")
+ self.assertEqual(pmb._create_nanoparticle_sites_positions(nanoparticle_tpl=tpl_zero), [])
+
+ def test_create_nanoparticle_zero_site_patch_branch(self):
+ """
+ Unit test: cover branch where a generated patch has zero sites.
+ """
+ pmb = self._build_pmb_with_particles()
+ pmb.define_nanoparticle(name="np",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+
+ original_helper = pmb._create_nanoparticle_sites_positions
+ pmb._create_nanoparticle_sites_positions = lambda nanoparticle_tpl: [
+ {"particle_name": "A", "positions": [], "number_of_sites": 0},
+ {"particle_name": "B", "positions": [[0.0, 0.0, 0.0]], "number_of_sites": 1},
+ ]
+ try:
+ created_ids = pmb.create_nanoparticle(name="np",
+ number_of_nanoparticles=1,
+ espresso_system=espresso_system,
+ list_core_particle_positions=[[1.0, 1.0, 1.0]],
+ fix=True)
+ finally:
+ pmb._create_nanoparticle_sites_positions = original_helper
+ self.assertEqual(created_ids, [0])
+
+ def test_nanoparticle_tools_standalone_functions(self):
+ """
+ Unit test: verify helper functions in ``nanoparticle_tools``.
+ """
+ points = nanoparticle_tools.uniform_distribution_sites_on_sphere(number_of_edges=1, tolerance=1e-6)
+ self.assertEqual(len(points), 1)
+ self.assertEqual(len(points[0]), 3)
+
+ distances = nanoparticle_tools.calculate_distance_vector_point([[0, 0, 0], [1, 0, 0]], [0, 0, 0])
+ self.assertEqual(distances[0], 0.0)
+ self.assertEqual(distances[1], 1.0)
+
+ _, patch = nanoparticle_tools.define_patch(points=[[0, 0, 0], [1, 0, 0], [0, 1, 0]],
+ central_point=[0, 0, 0],
+ patch_size=2)
+ self.assertEqual(len(patch), 2)
+
+ self.assertEqual(nanoparticle_tools.check_patch_overlaps(sites_positions=[[(0, 0, 0)], [(1, 0, 0)]],
+ number_patches=2), 0)
+ with self.assertRaises(ValueError):
+ nanoparticle_tools.check_patch_overlaps(sites_positions=[[(0, 0, 0)], [(0, 0, 0)]],
+ number_patches=2)
+
+ mean_d, std_d, err_d = nanoparticle_tools.calculate_distance_between_points_on_sphere(
+ points=[[(1.0, 0.0, 0.0),
+ (0.0, 1.0, 0.0),
+ (0.0, 0.0, 1.0),
+ (-1.0, 0.0, 0.0),
+ (0.0, -1.0, 0.0),
+ (0.0, 0.0, -1.0),
+ (0.707, 0.707, 0.0)]]
+ )
+ self.assertGreaterEqual(mean_d, 0.0)
+ self.assertGreaterEqual(std_d, 0.0)
+ self.assertGreaterEqual(err_d, 0.0)
+
+ dipole_vec, dipole_mag = nanoparticle_tools.calculate_dipole_moment(charges=[1.0, -1.0],
+ positions=[[1.0, 0.0, 0.0],
+ [0.0, 1.0, 0.0]])
+ self.assertEqual(len(dipole_vec), 3)
+ self.assertGreaterEqual(dipole_mag, 0.0)
+
+ q_tensor, q_mag, eigenvalues = nanoparticle_tools.calculate_quadrupole_moment(
+ charges=[1.0, -1.0],
+ positions=[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]],
+ )
+ self.assertEqual(q_tensor.shape, (3, 3))
+ self.assertEqual(len(eigenvalues), 3)
+ self.assertGreaterEqual(q_mag, 0.0)
+
+ def test_nanoparticle_template_edge_cases_and_manager_paths(self):
+ """
+ Unit test: verify nanoparticle-template edge cases and manager paths.
+
+ Notes:
+ - Covers template validation errors and nanoparticle-specific manager
+ collection logic.
+ """
+ pmb = self._build_pmb_with_particles()
+ pmb.define_nanoparticle(name="np",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+
+ # Manager nanoparticle branch in _collect_particle_templates and _get_templates_df
+ particle_counts = pmb.db.get_particle_templates_under(template_name="np",
+ pmb_type="nanoparticle",
+ return_counts=True)
+ self.assertEqual(set(particle_counts.keys()), {"core", "A", "B"})
+ self.assertFalse(pmb.get_templates_df("nanoparticle").empty)
+
+ # NanoparticleTemplate validations (fraction and number of patches)
+ tpl = pmb.db.get_template(pmb_type="nanoparticle", name="np")
+ with self.assertRaises(ValueError):
+ tpl.copy(update={"fraction_primary_sites": -0.1}).calculate_nanoparticle_properties(pmb)
+ with self.assertRaises(ValueError):
+ tpl.copy(update={"number_of_patches_of_primary_sites": 0}).calculate_nanoparticle_properties(pmb)
+
+ # Cover state_name is None path for _get_initial_state_charge_number
+ pmb.db._register_template(ParticleTemplate(name="A_no_init",
+ sigma=PintQuantity.from_quantity(1.0 * pmb.units.reduced_length, "length", pmb.units),
+ epsilon=PintQuantity.from_quantity(1.0 * pmb.units.reduced_energy, "energy", pmb.units),
+ cutoff=PintQuantity.from_quantity(1.2 * pmb.units.reduced_length, "length", pmb.units),
+ offset=PintQuantity.from_quantity(0.0 * pmb.units.reduced_length, "length", pmb.units),
+ initial_state=None))
+ pmb.define_particle_states(particle_name="A_no_init",
+ states=[{"name": "A_no_init_state", "z": 0}])
+ pmb.define_nanoparticle(name="np_no_init_site",
+ core_particle_name="core",
+ total_number_of_sites=5,
+ primary_site_particle_name="A_no_init",
+ fraction_primary_sites=1.0,
+ number_of_patches_of_primary_sites=1,
+ secondary_site_particle_name=None)
+ tpl_no_init = pmb.db.get_template(pmb_type="nanoparticle", name="np_no_init_site")
+ self.assertIn("total_number_of_sites", tpl_no_init.calculate_nanoparticle_properties(pmb))
+
+ # Cover core initial_state is None path
+ pmb.db._register_template(ParticleTemplate(name="core_no_init",
+ sigma=PintQuantity.from_quantity(4.0 * pmb.units.reduced_length, "length", pmb.units),
+ epsilon=PintQuantity.from_quantity(1.0 * pmb.units.reduced_energy, "energy", pmb.units),
+ cutoff=PintQuantity.from_quantity(4.2 * pmb.units.reduced_length, "length", pmb.units),
+ offset=PintQuantity.from_quantity(0.0 * pmb.units.reduced_length, "length", pmb.units),
+ initial_state=None))
+ pmb.define_nanoparticle(name="np_bad_core",
+ core_particle_name="core_no_init",
+ total_number_of_sites=5,
+ primary_site_particle_name="A",
+ fraction_primary_sites=1.0,
+ number_of_patches_of_primary_sites=1,
+ secondary_site_particle_name=None)
+ tpl_bad_core = pmb.db.get_template(pmb_type="nanoparticle", name="np_bad_core")
+ with self.assertRaises(ValueError):
+ tpl_bad_core.calculate_nanoparticle_properties(pmb)
+
+ # Cover no-particle-state branch for _get_initial_state_charge_number
+ pmb.db._register_template(ParticleTemplate(name="A_no_states",
+ sigma=PintQuantity.from_quantity(1.0 * pmb.units.reduced_length, "length", pmb.units),
+ epsilon=PintQuantity.from_quantity(1.0 * pmb.units.reduced_energy, "energy", pmb.units),
+ cutoff=PintQuantity.from_quantity(1.2 * pmb.units.reduced_length, "length", pmb.units),
+ offset=PintQuantity.from_quantity(0.0 * pmb.units.reduced_length, "length", pmb.units),
+ initial_state=None))
+ pmb.define_nanoparticle(name="np_no_states_site",
+ core_particle_name="core",
+ total_number_of_sites=5,
+ primary_site_particle_name="A_no_states",
+ fraction_primary_sites=1.0,
+ number_of_patches_of_primary_sites=1,
+ secondary_site_particle_name=None)
+ tpl_no_states = pmb.db.get_template(pmb_type="nanoparticle", name="np_no_states_site")
+ with self.assertRaises(ValueError):
+ tpl_no_states.calculate_nanoparticle_properties(pmb)
+
+ def test_nanoparticle_instance_and_io_roundtrip(self):
+ """
+ Unit test: verify nanoparticle instance validation and I/O roundtrip.
+
+ Notes:
+ - Confirms serialized nanoparticle templates and instances load back
+ correctly from disk.
+ """
+ pmb = self._build_pmb_with_particles()
+ pmb.define_nanoparticle(name="np",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+
+ with self.assertRaises(ValueError):
+ NanoparticleInstance(name="np", molecule_id=-1)
+
+ pmb.create_nanoparticle(name="np",
+ number_of_nanoparticles=1,
+ espresso_system=espresso_system,
+ fix=True)
+
+ new_pmb = pyMBE.pymbe_library(seed=43)
+ with tempfile.TemporaryDirectory() as tmp_directory:
+ pmb.save_database(tmp_directory)
+ new_pmb.load_database(tmp_directory)
+
+ self.assertFalse(new_pmb.get_templates_df(pmb_type="nanoparticle").empty)
+ self.assertFalse(new_pmb.get_instances_df(pmb_type="nanoparticle").empty)
+
+ def test_get_nanoparticle_properties(self):
+ """
+ Unit test: verify ``get_nanoparticle_properties`` returns the correct properties.
+
+ Notes:
+ - Checks that the function returns the same result as calling
+ ``calculate_nanoparticle_properties`` directly on the template.
+ - Verifies that all expected keys are present in the returned dict.
+ - Verifies that a non-existent name raises an error.
+ """
+ pmb = self._build_pmb_with_particles()
+ pmb.define_nanoparticle(name="np",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+
+ properties = nanoparticle_tools.get_nanoparticle_properties(pmb, "np")
+
+ expected_keys = ["nanoparticle_surface_area",
+ "nanoparticle_volume",
+ "total_number_of_sites",
+ "real_surface_density_of_sites",
+ "number_of_primary_sites",
+ "number_of_primary_sites_per_patch",
+ "number_of_secondary_sites",
+ "real_fraction_primary_sites",
+ "primary_site_charge_number",
+ "secondary_site_charge_number",
+ "total_charge",
+ "surface_charge_density",
+ "volume_charge_density"]
+ for key in expected_keys:
+ self.assertIn(key, properties)
+
+ self.assertEqual(properties["total_number_of_sites"], 10)
+
+ tpl = pmb.db.get_template(pmb_type="nanoparticle", name="np")
+ direct = tpl.calculate_nanoparticle_properties(pmb)
+ self.assertEqual(properties, direct)
+
+ with self.assertRaises(Exception):
+ nanoparticle_tools.get_nanoparticle_properties(pmb, "nonexistent_np")
+
+ def test_print_nanoparticle_properties(self):
+ """
+ Unit test: verify ``print_nanoparticle_properties`` prints all keys.
+
+ Notes:
+ - Checks that every key in the properties dict appears in stdout.
+ - Verifies that the nanoparticle name appears in the header.
+ - Verifies that the function works without a name argument.
+ """
+ pmb = self._build_pmb_with_particles()
+ pmb.define_nanoparticle(name="np",
+ core_particle_name="core",
+ total_number_of_sites=10,
+ primary_site_particle_name="A",
+ fraction_primary_sites=0.5,
+ number_of_patches_of_primary_sites=2,
+ secondary_site_particle_name="B")
+
+ properties = nanoparticle_tools.get_nanoparticle_properties(pmb, "np")
+
+ import io
+ import sys
+ captured = io.StringIO()
+ sys.stdout = captured
+ try:
+ nanoparticle_tools.print_nanoparticle_properties(properties, name="np")
+ finally:
+ sys.stdout = sys.__stdout__
+ output = captured.getvalue()
+
+ self.assertIn("np", output)
+ for key in properties:
+ self.assertIn(key, output)
+
+ # Also verify it runs without a name (no error, generic header)
+ captured2 = io.StringIO()
+ sys.stdout = captured2
+ try:
+ nanoparticle_tools.print_nanoparticle_properties(properties)
+ finally:
+ sys.stdout = sys.__stdout__
+ self.assertGreater(len(captured2.getvalue()), 0)
+
+
+if __name__ == "__main__":
+ ut.main()
diff --git a/testsuite/test_handy_functions.py b/testsuite/test_handy_functions.py
index 61ad4fce..8c343531 100644
--- a/testsuite/test_handy_functions.py
+++ b/testsuite/test_handy_functions.py
@@ -27,6 +27,7 @@
import io
import re
import numpy as np
+from unittest.mock import patch
# Create an in-memory log stream
log_stream = io.StringIO()
logging.basicConfig(level=logging.INFO,
@@ -307,5 +308,62 @@ def test_setup_electrostatics(self):
second=electrostatics_inputs["params"]["r_cut"],
msg="lib.handy_functions.setup_electrostatic_interactions sets up the wrong cut-off for the DH method")
+ def test_generate_lattice_positions(self):
+ """Test :func:`lib.handy_functions.generate_lattice_positions`."""
+ # Basic fcc generation in a box
+ positions_fcc = hf.generate_lattice_positions(lattice_type="fcc",
+ number_of_sites=10,
+ box_length=10.0)
+ self.assertEqual(len(positions_fcc), 10)
+ for pos in positions_fcc:
+ self.assertEqual(len(pos), 3)
+ self.assertTrue(all(0.0 <= value <= 10.0 for value in pos))
+
+ # bcc generation with explicit lattice constant and origin
+ origin = [1.0, 2.0, 3.0]
+ positions_bcc = hf.generate_lattice_positions(lattice_type="bcc",
+ number_of_sites=3,
+ lattice_constant=2.0,
+ origin=origin)
+ self.assertEqual(len(positions_bcc), 3)
+ self.assertListEqual(positions_bcc[0], origin)
+
+ # simple cubic and empty request
+ positions_sc = hf.generate_lattice_positions(lattice_type="sc",
+ number_of_sites=2,
+ lattice_constant=1.5)
+ self.assertEqual(len(positions_sc), 2)
+ self.assertEqual(hf.generate_lattice_positions(lattice_type="sc",
+ number_of_sites=0), [])
+
+ def test_generate_lattice_positions_exceptions(self):
+ """Test exceptions in :func:`lib.handy_functions.generate_lattice_positions`."""
+ with self.assertRaises(ValueError):
+ hf.generate_lattice_positions(lattice_type="hcp", number_of_sites=4)
+ with self.assertRaises(ValueError):
+ hf.generate_lattice_positions(lattice_type="fcc", number_of_sites=-1)
+ with self.assertRaises(ValueError):
+ hf.generate_lattice_positions(lattice_type="fcc", number_of_sites=1, box_length=0.0)
+ with self.assertRaises(ValueError):
+ hf.generate_lattice_positions(lattice_type="fcc", number_of_sites=1, lattice_constant=0.0)
+ with self.assertRaises(ValueError):
+ hf.generate_lattice_positions(lattice_type="fcc", number_of_sites=1, origin=[0.0, 1.0])
+
+ def test_generate_lattice_positions_internal_branches(self):
+ """Test internal branches in :func:`lib.handy_functions.generate_lattice_positions`."""
+ # Force n_cells <= 0 branch
+ with patch("pyMBE.lib.handy_functions.np.ceil", return_value=0):
+ positions = hf.generate_lattice_positions(lattice_type="sc",
+ number_of_sites=1,
+ lattice_constant=1.0)
+ self.assertEqual(len(positions), 1)
+
+ # Force fallback return path after nested loops
+ with patch("builtins.range", side_effect=lambda *args: []):
+ positions = hf.generate_lattice_positions(lattice_type="fcc",
+ number_of_sites=2,
+ lattice_constant=1.0)
+ self.assertEqual(positions, [])
+
if __name__ == "__main__":
- ut.main()
\ No newline at end of file
+ ut.main()