From 2786c10082ddf30a13f635ba83a53ddf88813f31 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Wed, 8 Apr 2026 17:43:40 +0200 Subject: [PATCH 01/13] Add support for angle potentials #143 --- pyMBE/pyMBE.py | 267 +++++++++++++++- pyMBE/storage/instances/angle.py | 58 ++++ pyMBE/storage/manager.py | 44 +++ pyMBE/storage/templates/angle.py | 109 +++++++ samples/build_angle.py | 118 +++++++ testsuite/angle_tests.py | 512 +++++++++++++++++++++++++++++++ 6 files changed, 1103 insertions(+), 5 deletions(-) create mode 100644 pyMBE/storage/instances/angle.py create mode 100644 pyMBE/storage/templates/angle.py create mode 100644 samples/build_angle.py create mode 100644 testsuite/angle_tests.py diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index f2bd213b..6134efc3 100644 --- a/pyMBE/pyMBE.py +++ b/pyMBE/pyMBE.py @@ -37,6 +37,7 @@ from pyMBE.storage.templates.protein import ProteinTemplate from pyMBE.storage.templates.hydrogel import HydrogelTemplate, HydrogelNode, HydrogelChain from pyMBE.storage.templates.bond import BondTemplate +from pyMBE.storage.templates.angle import AngleTemplate from pyMBE.storage.templates.lj import LJInteractionTemplate ## Instances from pyMBE.storage.instances.particle import ParticleInstance @@ -45,6 +46,7 @@ from pyMBE.storage.instances.peptide import PeptideInstance from pyMBE.storage.instances.protein import ProteinInstance from pyMBE.storage.instances.bond import BondInstance +from pyMBE.storage.instances.angle import AngleInstance from pyMBE.storage.instances.hydrogel import HydrogelInstance ## Reactions from pyMBE.storage.reactions.reaction import Reaction, ReactionParticipant @@ -956,7 +958,7 @@ def create_hydrogel(self, name, espresso_system, use_default_bond=False): assembly_id=assembly_id)) return assembly_id - def create_molecule(self, name, number_of_molecules, espresso_system, list_of_first_residue_positions=None, backbone_vector=None, use_default_bond=False, reverse_residue_order = False): + def create_molecule(self, name, number_of_molecules, espresso_system, list_of_first_residue_positions=None, backbone_vector=None, use_default_bond=False, reverse_residue_order = False, gen_angle=False): """ Creates instances of a given molecule template name into ESPResSo. @@ -1095,6 +1097,11 @@ def create_molecule(self, name, number_of_molecules, espresso_system, list_of_fi inst = PeptideInstance(name=name, molecule_id=molecule_id) self.db._register_instance(inst) + if gen_angle: + self._generate_angles_for_entity( + espresso_system=espresso_system, + entity_id=molecule_id, + entity_id_col='molecule_id') first_residue = True pos_index+=1 molecule_ids.append(molecule_id) @@ -1246,7 +1253,7 @@ def create_protein(self, name, number_of_proteins, espresso_system, topology_dic mol_ids.append(molecule_id) return mol_ids - def create_residue(self, name, espresso_system, central_bead_position=None,use_default_bond=False, backbone_vector=None): + def create_residue(self, name, espresso_system, central_bead_position=None,use_default_bond=False, backbone_vector=None, gen_angle=False): """ Creates a residue into ESPResSo. @@ -1370,8 +1377,12 @@ def create_residue(self, name, espresso_system, central_bead_position=None,use_d self.create_bond(particle_id1=central_bead_id, particle_id2=side_chain_beads_ids[0], espresso_system=espresso_system, - use_default_bond=use_default_bond) - return residue_id + use_default_bond=use_default_bond) + if gen_angle: + self._generate_angles_for_entity(espresso_system=espresso_system, + entity_id=residue_id, + entity_id_col="residue_id") + return residue_id def define_bond(self, bond_type, bond_parameters, particle_pairs): """ @@ -1452,7 +1463,253 @@ def define_default_bond(self, bond_type, bond_parameters): bond_type=bond_type) tpl.name = "default" self.db._register_template(tpl) - + + def define_angle(self, angle_type, angle_parameters, particle_triplets): + """ + Defines angle potential templates for each particle triplet in `particle_triplets`. + + Args: + angle_type ('str'): + Type of angle potential. Supported: "harmonic", "cosine", "harmonic_cosine". + + angle_parameters ('dict'): + Parameters of the angle potential. Must contain: + - "k" ('pint.Quantity'): Bending stiffness with dimensions of energy. + - "phi_0" ('float'): Equilibrium angle in radians. + + particle_triplets ('list[tuple[str,str,str]]'): + List of (side_particle1, central_particle, side_particle2) triplets. + """ + valid_angle_types = ["harmonic", "cosine", "harmonic_cosine"] + if angle_type not in valid_angle_types: + raise NotImplementedError(f"Angle potential type '{angle_type}' currently not implemented in pyMBE, accepted types are {valid_angle_types}") + + if "k" not in angle_parameters: + raise ValueError("Magnitude of the angle potential (k) is missing") + if "phi_0" not in angle_parameters: + raise ValueError("Equilibrium angle (phi_0) is missing") + + parameters_tpl = { + "k": PintQuantity.from_quantity(q=angle_parameters["k"], + expected_dimension="energy", + ureg=self.units), + "phi_0": PintQuantity.from_quantity(q=angle_parameters["phi_0"], + expected_dimension="dimensionless", + ureg=self.units), + } + + angle_names = [] + for side1, central, side2 in particle_triplets: + tpl = AngleTemplate(side_particle1=side1, + central_particle=central, + side_particle2=side2, + parameters=parameters_tpl, + angle_type=angle_type) + tpl._make_name() + if tpl.name in angle_names: + raise RuntimeError(f"Angle {tpl.name} has already been defined, please check the list of particle triplets") + angle_names.append(tpl.name) + self.db._register_template(tpl) + + def define_default_angle(self, angle_type, angle_parameters): + """ + Defines an angle template as a "default" template in the pyMBE database. + + Args: + angle_type ('str'): + Type of angle potential. Supported: "harmonic", "cosine", "harmonic_cosine". + + angle_parameters ('dict'): + Parameters of the angle potential (k, phi_0). + """ + valid_angle_types = ["harmonic", "cosine", "harmonic_cosine"] + if angle_type not in valid_angle_types: + raise NotImplementedError(f"Angle potential type '{angle_type}' currently not implemented in pyMBE, accepted types are {valid_angle_types}") + + if "k" not in angle_parameters: + raise ValueError("Magnitude of the angle potential (k) is missing") + if "phi_0" not in angle_parameters: + raise ValueError("Equilibrium angle (phi_0) is missing") + + parameters_tpl = { + "k": PintQuantity.from_quantity(q=angle_parameters["k"], + expected_dimension="energy", + ureg=self.units), + "phi_0": PintQuantity.from_quantity(q=angle_parameters["phi_0"], + expected_dimension="dimensionless", + ureg=self.units), + } + tpl = AngleTemplate(parameters=parameters_tpl, + angle_type=angle_type) + tpl.name = "default" + self.db._register_template(tpl) + + def create_angle(self, particle_id1, particle_id2, particle_id3, espresso_system, use_default_angle=False): + """ + Creates an angle between three particle instances in an ESPResSo system + and registers it in the pyMBE database. + + Args: + particle_id1 ('int'): ID of the first side particle. + particle_id2 ('int'): ID of the central particle. + particle_id3 ('int'): ID of the second side particle. + espresso_system ('espressomd.system.System'): ESPResSo system. + use_default_angle ('bool', optional): If True, use the default angle if no specific one is found. + """ + particle_inst_1 = self.db.get_instance(pmb_type="particle", instance_id=particle_id1) + particle_inst_2 = self.db.get_instance(pmb_type="particle", instance_id=particle_id2) + particle_inst_3 = self.db.get_instance(pmb_type="particle", instance_id=particle_id3) + + # Verify that bonds exist between side particles and central particle + bond_instances = self.db.get_instances(pmb_type="bond") + bonded_pairs = set() + for bond in bond_instances.values(): + pair = frozenset([bond.particle_id1, bond.particle_id2]) + bonded_pairs.add(pair) + if frozenset([particle_id1, particle_id2]) not in bonded_pairs: + raise ValueError(f"Cannot create angle: no bond exists between particle {particle_id1} and central particle {particle_id2}.") + if frozenset([particle_id3, particle_id2]) not in bonded_pairs: + raise ValueError(f"Cannot create angle: no bond exists between particle {particle_id3} and central particle {particle_id2}.") + + angle_tpl = self.get_angle_template(side_name1=particle_inst_1.name, + central_name=particle_inst_2.name, + side_name2=particle_inst_3.name, + use_default_angle=use_default_angle) + angle_inst = self._get_espresso_angle_instance(angle_template=angle_tpl, espresso_system=espresso_system) + + # ESPResSo angle bonds are added to the central particle + espresso_system.part.by_id(particle_id2).add_bond((angle_inst, particle_id1, particle_id3)) + + angle_id = self.db._propose_instance_id(pmb_type="angle") + pmb_angle_instance = AngleInstance(angle_id=angle_id, + name=angle_tpl.name, + particle_id1=particle_id1, + particle_id2=particle_id2, + particle_id3=particle_id3) + self.db._register_instance(instance=pmb_angle_instance) + + def get_angle_template(self, side_name1, central_name, side_name2, use_default_angle=False): + """ + Retrieves an angle template connecting three particle templates. + + Args: + side_name1 ('str'): Name of the first side particle. + central_name ('str'): Name of the central particle. + side_name2 ('str'): Name of the second side particle. + use_default_angle ('bool', optional): If True, fall back to the default angle template. + + Returns: + ('AngleTemplate'): The matching angle template. + """ + angle_key = AngleTemplate.make_angle_key(side1=side_name1, central=central_name, side2=side_name2) + try: + return self.db.get_template(name=angle_key, pmb_type="angle") + except ValueError: + pass + + if use_default_angle: + return self.db.get_template(name="default", pmb_type="angle") + + raise ValueError(f"No angle template found for '{side_name1}-{central_name}-{side_name2}', and default angles are deactivated.") + + def _get_espresso_angle_instance(self, angle_template, espresso_system): + """ + Retrieve or create an angle interaction in an ESPResSo system for a given angle template. + + Args: + angle_template ('AngleTemplate'): The angle template to use. + espresso_system ('espressomd.system.System'): ESPResSo system. + + Returns: + ('espressomd.interactions.BondedInteraction'): The ESPResSo angle interaction object. + """ + if angle_template.name in self.db.espresso_angle_instances: + return self.db.espresso_angle_instances[angle_template.name] + angle_inst = self._create_espresso_angle_instance(angle_type=angle_template.angle_type, + angle_parameters=angle_template.get_parameters(self.units)) + self.db.espresso_angle_instances[angle_template.name] = angle_inst + espresso_system.bonded_inter.add(angle_inst) + return angle_inst + + def _create_espresso_angle_instance(self, angle_type, angle_parameters): + """ + Creates an ESPResSo angle interaction object. + + Args: + angle_type ('str'): Type of angle potential ("harmonic", "cosine", "harmonic_cosine"). + angle_parameters ('dict'): Parameters of the angle potential (k, phi_0). + + Returns: + ('espressomd.interactions.BondedInteraction'): The ESPResSo angle interaction object. + """ + from espressomd import interactions + + k = angle_parameters["k"].m_as("reduced_energy") + phi_0 = float(angle_parameters["phi_0"].magnitude) + + if angle_type == "harmonic": + return interactions.AngleHarmonic(bend=k, phi0=phi_0) + elif angle_type == "cosine": + return interactions.AngleCosine(bend=k, phi0=phi_0) + elif angle_type == "harmonic_cosine": + return interactions.AngleCossquare(bend=k, phi0=phi_0) + else: + raise NotImplementedError(f"Angle potential type '{angle_type}' currently not implemented in pyMBE") + + def _generate_angles_for_entity(self, espresso_system, entity_id, entity_id_col): + """ + Auto-generates angles from bond topology for an entity (molecule or residue). + + For each particle in the entity that has two or more bonded neighbors, + this method finds all neighbor pairs and applies any matching angle potential. + + Args: + espresso_system ('espressomd.system.System'): ESPResSo system. + entity_id ('int'): The molecule_id or residue_id to generate angles for. + entity_id_col ('str'): Either "molecule_id" or "residue_id". + """ + # Get all particle IDs for this entity + particle_ids = self.db._find_instance_ids_by_attribute(pmb_type="particle", + attribute=entity_id_col, + value=entity_id) + if not particle_ids: + return + + # Build neighbor map from bond instances + neighbors = {pid: set() for pid in particle_ids} + pid_set = set(particle_ids) + bond_instances = self.db.get_instances(pmb_type="bond") + for bond in bond_instances.values(): + i, j = bond.particle_id1, bond.particle_id2 + if i in pid_set and j in pid_set: + neighbors[i].add(j) + neighbors[j].add(i) + + # For each particle with 2+ neighbors, generate angles + for j in particle_ids: + nbs = sorted(neighbors[j]) + if len(nbs) < 2: + continue + + central_name = self.db.get_instance(pmb_type="particle", instance_id=j).name + + for idx_i in range(len(nbs)): + for idx_k in range(idx_i + 1, len(nbs)): + i = nbs[idx_i] + k = nbs[idx_k] + side0_name = self.db.get_instance(pmb_type="particle", instance_id=i).name + side1_name = self.db.get_instance(pmb_type="particle", instance_id=k).name + + try: + self.create_angle(particle_id1=i, + particle_id2=j, + particle_id3=k, + espresso_system=espresso_system, + use_default_angle=True) + except ValueError: + # No angle template defined for this triplet — skip + continue + def define_hydrogel(self, name, node_map, chain_map): """ Defines a hydrogel template in the pyMBE database. diff --git a/pyMBE/storage/instances/angle.py b/pyMBE/storage/instances/angle.py new file mode 100644 index 00000000..34adbb80 --- /dev/null +++ b/pyMBE/storage/instances/angle.py @@ -0,0 +1,58 @@ +# +# 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 typing import Literal +from pyMBE.storage.base_type import PMBBaseModel +from pydantic import validator + +class AngleInstance(PMBBaseModel): + """ + Instance representation of an angle between three particles. + + Attributes: + pmb_type ('Literal["angle"]'): + Fixed identifier set to ``"angle"`` for all angle instances. + + angle_id ('int'): + Unique non-negative integer identifying this angle instance. + + name ('str'): + Name of the angle template from which this instance was created. + + particle_id1 ('int'): + ID of the first side particle. + + particle_id2 ('int'): + ID of the central particle. + + particle_id3 ('int'): + ID of the second side particle. + """ + pmb_type: Literal["angle"] = "angle" + angle_id: int + name: str + particle_id1: int + particle_id2: int + particle_id3: int + + @validator("angle_id", "particle_id1", "particle_id2", "particle_id3") + def validate_non_negative_int(cls, value, field): + if value < 0: + raise ValueError(f"{field.name} must be a non-negative integer.") + return value diff --git a/pyMBE/storage/manager.py b/pyMBE/storage/manager.py index f4b3f404..8131303e 100644 --- a/pyMBE/storage/manager.py +++ b/pyMBE/storage/manager.py @@ -24,10 +24,12 @@ from pyMBE.storage.templates.residue import ResidueTemplate from pyMBE.storage.templates.molecule import MoleculeTemplate from pyMBE.storage.templates.bond import BondTemplate +from pyMBE.storage.templates.angle import AngleTemplate from pyMBE.storage.instances.particle import ParticleInstance from pyMBE.storage.instances.residue import ResidueInstance from pyMBE.storage.instances.molecule import MoleculeInstance from pyMBE.storage.instances.bond import BondInstance +from pyMBE.storage.instances.angle import AngleInstance from pyMBE.storage.reactions.reaction import Reaction from pyMBE.storage.templates.peptide import PeptideTemplate from pyMBE.storage.instances.peptide import PeptideInstance @@ -90,6 +92,7 @@ def __init__(self,units): self._assembly_like_types = ["hydrogel"] self._pmb_types = ["particle", "residue"] + self._molecule_like_types + self._assembly_like_types self.espresso_bond_instances= {} + self.espresso_angle_instances= {} def _collect_particle_templates(self, name, pmb_type): """ @@ -176,6 +179,28 @@ def _delete_bonds_of_particle(self, pid): if "bond" in self._instances and not self._instances["bond"]: del self._instances["bond"] + def _delete_angles_of_particle(self, pid): + """ + Delete all angle instances involving a given particle instance. + + Args: + pid ('int'): + The particle ID whose associated angles should be deleted. + + Notes: + - If no `"angle"` instances are present in the database, the method + exits immediately. + - This method does not raise errors if no angles involve the particle. + - It is intended for internal use by cascade-deletion routines. + """ + if "angle" not in self._instances: + return + angles_to_delete = [a_id for a_id, a in list(self._instances["angle"].items()) if a.particle_id1 == pid or a.particle_id2 == pid or a.particle_id3 == pid] + for a_id in angles_to_delete: + del self._instances["angle"][a_id] + if "angle" in self._instances and not self._instances["angle"]: + del self._instances["angle"] + def _find_instance_ids_by_attribute(self, pmb_type, attribute, value): """ Return a list of instance IDs for a given pmb_type where a given attribute @@ -355,6 +380,17 @@ def _get_templates_df(self, pmb_type): "particle_name1": tpl.particle_name1, "particle_name2": tpl.particle_name2, "parameters": parameters}) + elif pmb_type == "angle": + parameters = {} + for key in tpl.parameters.keys(): + parameters[key] = tpl.parameters[key].to_quantity(self._units) + rows.append({"pmb_type": tpl.pmb_type, + "name": tpl.name, + "angle_type": tpl.angle_type, + "side_particle1": tpl.side_particle1, + "central_particle": tpl.central_particle, + "side_particle2": tpl.side_particle2, + "parameters": parameters}) else: # Generic representation for other types rows.append(tpl.dict()) @@ -428,6 +464,9 @@ def _register_instance(self, instance): elif isinstance(instance, BondInstance): pmb_type = "bond" iid = instance.bond_id + elif isinstance(instance, AngleInstance): + pmb_type = "angle" + iid = instance.angle_id elif isinstance(instance, HydrogelInstance): pmb_type = "hydrogel" iid = instance.assembly_id @@ -679,6 +718,7 @@ def delete_instance(self, pmb_type, instance_id): # --- Delete children of PARTICLE (only bonds) --- if pmb_type == "particle": self._delete_bonds_of_particle(instance_id) + self._delete_angles_of_particle(instance_id) # =============== FINAL DELETION STEP ====================== del self._instances[pmb_type][instance_id] if not self._instances[pmb_type]: @@ -754,6 +794,10 @@ def delete_template(self, pmb_type, name): if pmb_type == "bond": if name in self.espresso_bond_instances.keys(): del self.espresso_bond_instances[name] + # if it is an angle template delete also stored espresso angle instances + if pmb_type == "angle": + if name in self.espresso_angle_instances.keys(): + del self.espresso_angle_instances[name] # Delete empty groups if not self._templates[pmb_type]: del self._templates[pmb_type] diff --git a/pyMBE/storage/templates/angle.py b/pyMBE/storage/templates/angle.py new file mode 100644 index 00000000..af5525d5 --- /dev/null +++ b/pyMBE/storage/templates/angle.py @@ -0,0 +1,109 @@ +# +# 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 typing import Dict, Literal, Optional +from ..base_type import PMBBaseModel +from ..pint_quantity import PintQuantity +from pydantic import Field + +class AngleTemplate(PMBBaseModel): + """ + Template defining an angular potential in the pyMBE database. + + Attributes: + pmb_type ('Literal["angle"]'): + Fixed type identifier for this template. Always "angle". + + name ('str'): + Unique name of the angle template, e.g., "A-B-C" where B is the central particle. + + angle_type ('str'): + Type of angle potential. Examples: "harmonic", "cosine", "harmonic_cosine". + + side_particle1 ('Optional[str]'): + Name of the first side particle in the angle triplet. + + central_particle ('Optional[str]'): + Name of the central particle in the angle triplet. + + side_particle2 ('Optional[str]'): + Name of the second side particle in the angle triplet. + + parameters ('Dict[str, PintQuantity]'): + Dictionary of angle parameters (e.g., k, phi_0). + + Notes: + - Values of the parameters are stored as PintQuantity objects for unit-aware calculations. + - The canonical name sorts the two side particles alphabetically while keeping + the central particle in the middle: ``"side_min-central-side_max"``. + """ + pmb_type: Literal["angle"] = "angle" + name: str = Field(default="default") + angle_type: str + side_particle1: Optional[str] = None + central_particle: Optional[str] = None + side_particle2: Optional[str] = None + parameters: Dict[str, PintQuantity] + + @classmethod + def make_angle_key(cls, side1, central, side2): + """Return a canonical name for an angle between three particle names. + + The two side particles are sorted alphabetically, with the central + particle kept in the middle position. + + Args: + side1 ('str'): + Name of the first side particle. + + central ('str'): + Name of the central particle. + + side2 ('str'): + Name of the second side particle. + + Returns: + ('str'): + Canonical angle name, e.g. "A-B-C". + """ + sides = sorted([side1, side2]) + return f"{sides[0]}-{central}-{sides[1]}" + + def _make_name(self): + """Create canonical name using particle names.""" + if not self.side_particle1 or not self.central_particle or not self.side_particle2: + raise RuntimeError("Cannot generate angle name: side_particle1, central_particle, or side_particle2 missing.") + self.name = self.make_angle_key(self.side_particle1, self.central_particle, self.side_particle2) + + def get_parameters(self, ureg): + """ + Retrieve the angle parameters as Pint `Quantity` objects. + + Args: + ureg ('pint.UnitRegistry'): + Pint unit registry used to reconstruct physical quantities from storage. + + Returns: + 'Dict[str, pint.Quantity]': + A dictionary mapping parameter names to their corresponding unit-aware Pint quantities. + """ + pint_parameters = {} + for parameter in self.parameters.keys(): + pint_parameters[parameter] = self.parameters[parameter].to_quantity(ureg) + return pint_parameters diff --git a/samples/build_angle.py b/samples/build_angle.py new file mode 100644 index 00000000..462eff1e --- /dev/null +++ b/samples/build_angle.py @@ -0,0 +1,118 @@ +# +# Example: demonstrating angle potential support in pyMBE. +# +# This script shows two complementary usages of `gen_angle=True`: +# 1. create_residue with gen_angle=True +# - builds a single residue with auto-generated angles +# 2. create_molecule with gen_angle=True +# - builds a chain of residues; angles are auto-generated across the +# full molecule, including angles that span residue boundaries +# +# Note: angle auto-generation is currently supported only in `create_residue` +# and `create_molecule`. Other higher-level builders (e.g. `create_hydrogel`) +# do not yet propagate `gen_angle`. +# + +import pyMBE +import numpy as np +import espressomd + +# ---------------------------------------------------------------------- +# Set up the ESPResSo system and pyMBE +# ---------------------------------------------------------------------- +espresso_system = espressomd.System(box_l=[20] * 3) +pmb = pyMBE.pymbe_library(seed=42) + +# ---------------------------------------------------------------------- +# Particle definitions +# ---------------------------------------------------------------------- +pmb.define_particle(name='A', + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + +pmb.define_particle(name='B', + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + +pmb.define_particle(name='C', + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + +# ---------------------------------------------------------------------- +# Default bond used everywhere +# ---------------------------------------------------------------------- +HARMONIC_bond_parameters = { + 'r_0': 0.4 * pmb.units.nm, + 'k': 40 * pmb.units('reduced_energy / reduced_length**2'), +} +pmb.define_default_bond(bond_type='harmonic', + bond_parameters=HARMONIC_bond_parameters) + +# ---------------------------------------------------------------------- +# Default angle used by both demos +# ---------------------------------------------------------------------- +default_angle_params = { + 'k': 50 * pmb.units('reduced_energy'), + 'phi_0': np.pi * pmb.units(''), +} +pmb.define_default_angle(angle_type="harmonic", + angle_parameters=default_angle_params) + + +def show_database(label): + print(f"\n############ {label} ############") + for pmb_type in ["particle", "residue", "molecule", "bond", "angle"]: + print(f"\n=== {pmb_type} templates ===") + print(pmb.get_templates_df(pmb_type=pmb_type)) + print(f"\n=== {pmb_type} instances ===") + print(pmb.get_instances_df(pmb_type=pmb_type)) + + +# ====================================================================== +# Demo 1: create_residue with gen_angle=True +# ====================================================================== +# Topology: central A bonded to side particles A, B, C. +# Auto-generated angles (central A has 3 neighbors): A-A-B, A-A-C, B-A-C +pmb.define_residue(name="Res_1", + central_bead="A", + side_chains=["A", "B", "C"]) + +residue_id = pmb.create_residue(name="Res_1", + espresso_system=espresso_system, + use_default_bond=True, + gen_angle=True) + +show_database("Demo 1: single residue with gen_angle=True") + +# ---------------------------------------------------------------------- +# Clean up: remove the residue (and its particles, bonds, angles) from +# both ESPResSo and the pyMBE database before running Demo 2. +# ---------------------------------------------------------------------- +pmb.delete_instances_in_system(instance_id=residue_id, + pmb_type="residue", + espresso_system=espresso_system) + +# ====================================================================== +# Demo 2: create_molecule with gen_angle=True +# ====================================================================== +# Topology: a chain of 4 Res_2 residues, each with central A and one +# side chain B. After auto-generation, angles include both intra-residue +# (A-A-B) and inter-residue (A-A-A) triplets. +pmb.define_residue(name="Res_2", + central_bead="A", + side_chains=["B"]) + +pmb.define_molecule(name="Mol_1", + residue_list=["Res_2"] * 4) + +pmb.create_molecule(name="Mol_1", + number_of_molecules=1, + espresso_system=espresso_system, + backbone_vector=[1.0, 0.0, 0.0], + use_default_bond=True, + gen_angle=True) + +show_database("Demo 2: linear molecule with gen_angle=True") diff --git a/testsuite/angle_tests.py b/testsuite/angle_tests.py new file mode 100644 index 00000000..8a986f88 --- /dev/null +++ b/testsuite/angle_tests.py @@ -0,0 +1,512 @@ +# +# Copyright (C) 2024-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 pyMBE and other libraries +import pyMBE +import numpy as np +import unittest as ut +import espressomd + +# Create an instance of the ESPResSo system +espresso_system = espressomd.System(box_l=[10]*3) + + +class Test(ut.TestCase): + def setUp(self): + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + + def define_templates(self, pmb): + pmb.define_particle(name='A', + z=0, + sigma=0.4*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + + pmb.define_particle(name='B', + z=0, + sigma=0.4*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + + pmb.define_particle(name='C', + z=0, + sigma=0.4*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + + self.harmonic_angle_params = { + 'k': 50 * pmb.units('reduced_energy'), + 'phi_0': np.pi * pmb.units(''), + } + + self.cosine_angle_params = { + 'k': 30 * pmb.units('reduced_energy'), + 'phi_0': np.pi / 2 * pmb.units(''), + } + + self.harmonic_bond_params = { + 'r_0': 0.4 * pmb.units.nm, + 'k': 400 * pmb.units('reduced_energy / reduced_length**2'), + } + + def get_angle_object(self, central_particle_id): + """ + Returns the angle interaction object stored in ESPResSo for a given central particle. + Angle bonds are added to the central particle in ESPResSo. + """ + bonds = espresso_system.part.by_id(central_particle_id).bonds + for bond_tuple in bonds: + bond_obj = bond_tuple[0] # Testing only one angle on the central particle + # Angle interactions involve 2 partner particles (the two sides) + if len(bond_tuple) == 3: + return bond_obj + return None + + def check_angle_setup(self, angle_object, input_parameters, angle_type): + """ + Checks that pyMBE sets up an angle interaction object correctly. + + Args: + angle_object: instance of an ESPResSo angle interaction object. + input_parameters('dict'): dictionary with the parameters for the angle potential. + angle_type('str'): label identifying the angle type. + """ + # Check that pyMBE stores the correct type of angle + type_map = { + "harmonic": "AngleHarmonic", + "cosine": "AngleCosine", + "harmonic_cosine": "AngleCossquare", + } + self.assertIn(type_map[angle_type].lower(), str(type(angle_object)).lower(), + msg=f"pyMBE does not store the correct type of angle interaction, expected {type_map[angle_type]}") + # Check that pyMBE defines the angle with the correct parameters + angle_params = angle_object.get_params() + np.testing.assert_almost_equal(actual=angle_params['bend'], + desired=input_parameters['k'].m_as('reduced_energy'), + verbose=True) + np.testing.assert_almost_equal(actual=angle_params['phi0'], + desired=float(input_parameters['phi_0'].magnitude), + verbose=True) + + def test_angle_setup_harmonic(self): + """ + Unit test to check the setup of harmonic angle potentials in pyMBE. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + # Define bonds between A-B and B-C + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B'], ['B', 'C']]) + + # Define harmonic angle for triplet A-B-C + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('A', 'B', 'C')]) + + # Create three particles: A, B, C + pid_A = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_C = pmb.create_particle(name="C", + espresso_system=espresso_system, + number_of_particles=1) + + # Create bonds: A-B and B-C + pmb.create_bond(particle_id1=pid_A[0], + particle_id2=pid_B[0], + espresso_system=espresso_system) + pmb.create_bond(particle_id1=pid_B[0], + particle_id2=pid_C[0], + espresso_system=espresso_system) + + # Create the angle: A-B-C (B is central) + pmb.create_angle(particle_id1=pid_A[0], + particle_id2=pid_B[0], + particle_id3=pid_C[0], + espresso_system=espresso_system) + + angle_object = self.get_angle_object(central_particle_id=pid_B[0]) + self.assertIsNotNone(angle_object, "No angle object found on the central particle") + + self.check_angle_setup(angle_object=angle_object, + input_parameters=self.harmonic_angle_params, + angle_type="harmonic") + + # Clean-up + for inst_id in pid_A + pid_B + pid_C: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + + def test_angle_setup_cosine(self): + """ + Unit test to check the setup of cosine angle potentials in pyMBE. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + # Define bonds between A-B + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B']]) + + # Define cosine angle for triplet A-B-A + pmb.define_angle(angle_type="cosine", + angle_parameters=self.cosine_angle_params, + particle_triplets=[('A', 'B', 'A')]) + + # Create three particles + pid_A1 = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_A2 = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + + # Create bonds: A1-B and B-A2 + pmb.create_bond(particle_id1=pid_A1[0], + particle_id2=pid_B[0], + espresso_system=espresso_system) + pmb.create_bond(particle_id1=pid_B[0], + particle_id2=pid_A2[0], + espresso_system=espresso_system) + + pmb.create_angle(particle_id1=pid_A1[0], + particle_id2=pid_B[0], + particle_id3=pid_A2[0], + espresso_system=espresso_system) + + angle_object = self.get_angle_object(central_particle_id=pid_B[0]) + self.assertIsNotNone(angle_object, "No angle object found on the central particle") + + self.check_angle_setup(angle_object=angle_object, + input_parameters=self.cosine_angle_params, + angle_type="cosine") + + # Clean-up + for inst_id in pid_A1 + pid_B + pid_A2: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + + def test_default_angle(self): + """ + Unit test to check the setup of default angle potentials. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + # Define bonds between A-B and B-C + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B'], ['B', 'C']]) + + # Define only a default angle (no specific triplet) + pmb.define_default_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params) + + # Create three particles + pid_A = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_C = pmb.create_particle(name="C", + espresso_system=espresso_system, + number_of_particles=1) + + # Create bonds: A-B and B-C + pmb.create_bond(particle_id1=pid_A[0], + particle_id2=pid_B[0], + espresso_system=espresso_system) + pmb.create_bond(particle_id1=pid_B[0], + particle_id2=pid_C[0], + espresso_system=espresso_system) + + # Create angle using default (no specific A-B-C template exists) + pmb.create_angle(particle_id1=pid_A[0], + particle_id2=pid_B[0], + particle_id3=pid_C[0], + espresso_system=espresso_system, + use_default_angle=True) + + angle_object = self.get_angle_object(central_particle_id=pid_B[0]) + self.assertIsNotNone(angle_object, "No angle object found on the central particle") + + self.check_angle_setup(angle_object=angle_object, + input_parameters=self.harmonic_angle_params, + angle_type="harmonic") + + # Clean-up + for inst_id in pid_A + pid_B + pid_C: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + + def test_specific_angle_over_default(self): + """ + Test that a specific angle template is preferred over the default. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + # Define bonds between A-B and B-C + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B'], ['B', 'C']]) + + # Define default angle with cosine parameters + pmb.define_default_angle(angle_type="cosine", + angle_parameters=self.cosine_angle_params) + + # Define specific angle for A-B-C with harmonic parameters + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('A', 'B', 'C')]) + + pid_A = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_C = pmb.create_particle(name="C", + espresso_system=espresso_system, + number_of_particles=1) + + # Create bonds: A-B and B-C + pmb.create_bond(particle_id1=pid_A[0], + particle_id2=pid_B[0], + espresso_system=espresso_system) + pmb.create_bond(particle_id1=pid_B[0], + particle_id2=pid_C[0], + espresso_system=espresso_system) + + # Should use the specific A-B-C template, not the default + pmb.create_angle(particle_id1=pid_A[0], + particle_id2=pid_B[0], + particle_id3=pid_C[0], + espresso_system=espresso_system, + use_default_angle=True) + + angle_object = self.get_angle_object(central_particle_id=pid_B[0]) + self.assertIsNotNone(angle_object) + + # Should be harmonic (specific), not cosine (default) + self.check_angle_setup(angle_object=angle_object, + input_parameters=self.harmonic_angle_params, + angle_type="harmonic") + + # Clean-up + for inst_id in pid_A + pid_B + pid_C: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + + def test_angle_raised_exceptions(self): + """ + Unit test to check that angle-related methods raise the correct exceptions. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + for callback in [pmb.define_angle, pmb.define_default_angle]: + with self.subTest(msg=f'using method {callback.__qualname__}()'): + self.check_angle_exceptions(callback, pmb) + + def check_angle_exceptions(self, callback, pmb): + # Check exception for unknown angle type + angle_type = 'quartic' + angle_params = {'k': 50 * pmb.units('reduced_energy'), + 'phi_0': np.pi * pmb.units('')} + + input_parameters = {"angle_type": angle_type, "angle_parameters": angle_params} + if callback == pmb.define_angle: + input_parameters["particle_triplets"] = [('A', 'B', 'C')] + + np.testing.assert_raises(NotImplementedError, callback, **input_parameters) + + # Check exception for missing k + angle_type = 'harmonic' + angle_params = {'phi_0': np.pi * pmb.units('')} + + input_parameters = {"angle_type": angle_type, "angle_parameters": angle_params} + if callback == pmb.define_angle: + input_parameters["particle_triplets"] = [('A', 'B', 'C')] + + np.testing.assert_raises(ValueError, callback, **input_parameters) + + # Check exception for missing phi_0 + angle_type = 'harmonic' + angle_params = {'k': 50 * pmb.units('reduced_energy')} + + input_parameters = {"angle_type": angle_type, "angle_parameters": angle_params} + if callback == pmb.define_angle: + input_parameters["particle_triplets"] = [('A', 'B', 'C')] + + np.testing.assert_raises(ValueError, callback, **input_parameters) + + # Check exception for duplicate triplet in define_angle + if callback == pmb.define_angle: + test = {"angle_type": "harmonic", + "angle_parameters": self.harmonic_angle_params, + "particle_triplets": [("A", "B", "A"), ("A", "B", "A")]} + + np.testing.assert_raises(RuntimeError, + pmb.define_angle, + **test) + + def test_sanity_get_angle_template(self): + """ + Tests that get_angle_template raises ValueError when no template is found. + """ + pmb = pyMBE.pymbe_library(seed=51) + inputs = {"side_name1": "X", + "central_name": "Y", + "side_name2": "Z"} + np.testing.assert_raises(ValueError, + pmb.get_angle_template, + **inputs) + + def test_canonical_angle_name(self): + """ + Test that angle names are canonical (side particles sorted alphabetically). + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + # Define angle with C-B-A order — should produce same template as A-B-C + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('C', 'B', 'A')]) + + # The canonical name should be A-B-C (sides sorted, central stays in middle) + tpl = pmb.get_angle_template(side_name1="A", + central_name="B", + side_name2="C") + self.assertEqual(tpl.name, "A-B-C") + + # Also retrievable with reversed side order + tpl2 = pmb.get_angle_template(side_name1="C", + central_name="B", + side_name2="A") + self.assertEqual(tpl2.name, "A-B-C") + + pmb.db.delete_templates(pmb_type="angle") + + def test_angle_without_bonds_raises_error(self): + """ + Test that creating an angle without bonds between the particles raises a ValueError. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + # Define angle template but no bond templates + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('A', 'B', 'C')]) + + # Create three particles without any bonds between them + pid_A = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_C = pmb.create_particle(name="C", + espresso_system=espresso_system, + number_of_particles=1) + + # Attempting to create angle without bonds should raise ValueError + with self.assertRaises(ValueError, msg="create_angle should raise ValueError when no bonds exist"): + pmb.create_angle(particle_id1=pid_A[0], + particle_id2=pid_B[0], + particle_id3=pid_C[0], + espresso_system=espresso_system) + + # Clean-up + for inst_id in pid_A + pid_B + pid_C: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + + def test_angle_with_partial_bonds_raises_error(self): + """ + Test that creating an angle with only one bond (missing the other) raises a ValueError. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + # Define bond only between A-B, not B-C + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B']]) + + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('A', 'B', 'C')]) + + pid_A = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_C = pmb.create_particle(name="C", + espresso_system=espresso_system, + number_of_particles=1) + + # Create only the A-B bond, not B-C + pmb.create_bond(particle_id1=pid_A[0], + particle_id2=pid_B[0], + espresso_system=espresso_system) + + # Should raise ValueError because B-C bond is missing + with self.assertRaises(ValueError, msg="create_angle should raise ValueError when a bond is missing"): + pmb.create_angle(particle_id1=pid_A[0], + particle_id2=pid_B[0], + particle_id3=pid_C[0], + espresso_system=espresso_system) + + # Clean-up + for inst_id in pid_A + pid_B + pid_C: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + + +if __name__ == '__main__': + ut.main() From 84e704e3a14eb9201b1fec9733c66024addec7bc Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Wed, 8 Apr 2026 18:02:11 +0200 Subject: [PATCH 02/13] unused variable removed --- pyMBE/pyMBE.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index 6134efc3..d915da78 100644 --- a/pyMBE/pyMBE.py +++ b/pyMBE/pyMBE.py @@ -1691,15 +1691,10 @@ def _generate_angles_for_entity(self, espresso_system, entity_id, entity_id_col) if len(nbs) < 2: continue - central_name = self.db.get_instance(pmb_type="particle", instance_id=j).name - for idx_i in range(len(nbs)): for idx_k in range(idx_i + 1, len(nbs)): i = nbs[idx_i] k = nbs[idx_k] - side0_name = self.db.get_instance(pmb_type="particle", instance_id=i).name - side1_name = self.db.get_instance(pmb_type="particle", instance_id=k).name - try: self.create_angle(particle_id1=i, particle_id2=j, From ffc61f00a6e6922264f406f99a53cd5786723341 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Tue, 21 Apr 2026 10:56:36 +0200 Subject: [PATCH 03/13] angle test coverage --- pyMBE/pyMBE.py | 2 - testsuite/CTestTestfile.cmake | 1 + testsuite/angle_tests.py | 278 ++++++++++++++++++++++++++++++++++ 3 files changed, 279 insertions(+), 2 deletions(-) diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index d915da78..f91dd244 100644 --- a/pyMBE/pyMBE.py +++ b/pyMBE/pyMBE.py @@ -1653,8 +1653,6 @@ def _create_espresso_angle_instance(self, angle_type, angle_parameters): return interactions.AngleCosine(bend=k, phi0=phi_0) elif angle_type == "harmonic_cosine": return interactions.AngleCossquare(bend=k, phi0=phi_0) - else: - raise NotImplementedError(f"Angle potential type '{angle_type}' currently not implemented in pyMBE") def _generate_angles_for_entity(self, espresso_system, entity_id, entity_id_col): """ diff --git a/testsuite/CTestTestfile.cmake b/testsuite/CTestTestfile.cmake index aa66fb39..4dfeea96 100644 --- a/testsuite/CTestTestfile.cmake +++ b/testsuite/CTestTestfile.cmake @@ -58,6 +58,7 @@ pymbe_add_test(PATH test_global_variables.py) pymbe_add_test(PATH lj_tests.py) pymbe_add_test(PATH set_particle_acidity_test.py) pymbe_add_test(PATH bond_tests.py) +pymbe_add_test(PATH angle_tests.py) pymbe_add_test(PATH generate_perpendicular_vectors_test.py) pymbe_add_test(PATH define_and_create_molecules_unit_tests.py) pymbe_add_test(PATH create_molecule_position_test.py) diff --git a/testsuite/angle_tests.py b/testsuite/angle_tests.py index 8a986f88..68a6afa8 100644 --- a/testsuite/angle_tests.py +++ b/testsuite/angle_tests.py @@ -57,6 +57,11 @@ def define_templates(self, pmb): 'phi_0': np.pi / 2 * pmb.units(''), } + self.harmonic_cosine_angle_params = { + 'k': 40 * pmb.units('reduced_energy'), + 'phi_0': np.pi / 3 * pmb.units(''), + } + self.harmonic_bond_params = { 'r_0': 0.4 * pmb.units.nm, 'k': 400 * pmb.units('reduced_energy / reduced_length**2'), @@ -214,6 +219,82 @@ def test_angle_setup_cosine(self): pmb.db.delete_templates(pmb_type="angle") pmb.db.delete_templates(pmb_type="bond") + def test_angle_setup_harmonic_cosine(self): + """ + Unit test to check the setup of harmonic-cosine angle potentials in pyMBE. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B'], ['B', 'C']]) + + pmb.define_angle(angle_type="harmonic_cosine", + angle_parameters=self.harmonic_cosine_angle_params, + particle_triplets=[('A', 'B', 'C')]) + + pid_A = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_C = pmb.create_particle(name="C", + espresso_system=espresso_system, + number_of_particles=1) + + pmb.create_bond(particle_id1=pid_A[0], + particle_id2=pid_B[0], + espresso_system=espresso_system) + pmb.create_bond(particle_id1=pid_B[0], + particle_id2=pid_C[0], + espresso_system=espresso_system) + + pmb.create_angle(particle_id1=pid_A[0], + particle_id2=pid_B[0], + particle_id3=pid_C[0], + espresso_system=espresso_system) + + angle_object = self.get_angle_object(central_particle_id=pid_B[0]) + self.assertIsNotNone(angle_object, "No angle object found on the central particle") + + self.check_angle_setup(angle_object=angle_object, + input_parameters=self.harmonic_cosine_angle_params, + angle_type="harmonic_cosine") + + for inst_id in pid_A + pid_B + pid_C: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + + def test_get_espresso_angle_instance_reuses_cached_object(self): + """ + Test that cached ESPResSo angle interactions are reused for the same template. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('A', 'B', 'C')]) + + angle_template = pmb.get_angle_template(side_name1="A", + central_name="B", + side_name2="C") + + first_angle_object = pmb._get_espresso_angle_instance(angle_template=angle_template, + espresso_system=espresso_system) + second_angle_object = pmb._get_espresso_angle_instance(angle_template=angle_template, + espresso_system=espresso_system) + + self.assertIs(first_angle_object, second_angle_object) + self.assertEqual(len(pmb.db.espresso_angle_instances), 1) + + pmb.db.delete_templates(pmb_type="angle") + def test_default_angle(self): """ Unit test to check the setup of default angle potentials. @@ -507,6 +588,203 @@ def test_angle_with_partial_bonds_raises_error(self): pmb.db.delete_templates(pmb_type="angle") pmb.db.delete_templates(pmb_type="bond") + def test_generate_angles_for_entity_without_particles_returns(self): + """ + Test that auto-generating angles for an empty entity returns without changes. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb._generate_angles_for_entity(espresso_system=espresso_system, + entity_id=999, + entity_id_col="residue_id") + + self.assertEqual(len(pmb.get_instances_df(pmb_type="angle")), 0) + + def test_generate_angles_for_entity_creates_angle_from_bonds(self): + """ + Test that auto-generating angles creates one angle from a bonded triplet. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B'], ['B', 'C']]) + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('A', 'B', 'C')]) + + pid_A = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_C = pmb.create_particle(name="C", + espresso_system=espresso_system, + number_of_particles=1) + + for particle_id in (pid_A[0], pid_B[0], pid_C[0]): + pmb.db._update_instance(instance_id=particle_id, + pmb_type="particle", + attribute="residue_id", + value=0) + + pmb.create_bond(particle_id1=pid_A[0], + particle_id2=pid_B[0], + espresso_system=espresso_system) + pmb.create_bond(particle_id1=pid_B[0], + particle_id2=pid_C[0], + espresso_system=espresso_system) + + pmb._generate_angles_for_entity(espresso_system=espresso_system, + entity_id=0, + entity_id_col="residue_id") + + angle_object = self.get_angle_object(central_particle_id=pid_B[0]) + self.assertIsNotNone(angle_object, "No angle object found on the central particle") + self.assertEqual(len(pmb.get_instances_df(pmb_type="angle")), 1) + + for inst_id in pid_A + pid_B + pid_C: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + + def test_generate_angles_for_entity_skips_missing_templates(self): + """ + Test that auto-generating angles skips bonded triplets without a matching angle template. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B'], ['B', 'C']]) + + pid_A = pmb.create_particle(name="A", + espresso_system=espresso_system, + number_of_particles=1) + pid_B = pmb.create_particle(name="B", + espresso_system=espresso_system, + number_of_particles=1) + pid_C = pmb.create_particle(name="C", + espresso_system=espresso_system, + number_of_particles=1) + + for particle_id in (pid_A[0], pid_B[0], pid_C[0]): + pmb.db._update_instance(instance_id=particle_id, + pmb_type="particle", + attribute="residue_id", + value=0) + + pmb.create_bond(particle_id1=pid_A[0], + particle_id2=pid_B[0], + espresso_system=espresso_system) + pmb.create_bond(particle_id1=pid_B[0], + particle_id2=pid_C[0], + espresso_system=espresso_system) + + pmb._generate_angles_for_entity(espresso_system=espresso_system, + entity_id=0, + entity_id_col="residue_id") + + self.assertIsNone(self.get_angle_object(central_particle_id=pid_B[0])) + self.assertEqual(len(pmb.get_instances_df(pmb_type="angle")), 0) + + for inst_id in pid_A + pid_B + pid_C: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="bond") + + def test_create_residue_with_gen_angle_generates_angles(self): + """ + Test that create_residue generates angles automatically when gen_angle is enabled. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B'], ['B', 'C']]) + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('A', 'B', 'C')]) + pmb.define_residue(name="R_angle", + central_bead="B", + side_chains=["A", "C"]) + + residue_id = pmb.create_residue(name="R_angle", + espresso_system=espresso_system, + gen_angle=True) + + particle_ids = pmb.db._find_instance_ids_by_attribute(pmb_type="particle", + attribute="residue_id", + value=residue_id) + central_particle_id = next(pid for pid in particle_ids + if pmb.db.get_instance(pmb_type="particle", + instance_id=pid).name == "B") + + angle_object = self.get_angle_object(central_particle_id=central_particle_id) + self.assertIsNotNone(angle_object, "No angle object found on the central particle") + self.assertEqual(len(pmb.get_instances_df(pmb_type="angle")), 1) + + pmb.delete_instances_in_system(instance_id=residue_id, + pmb_type="residue", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + + def test_create_molecule_with_gen_angle_generates_angles(self): + """ + Test that create_molecule generates backbone angles automatically when gen_angle is enabled. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_bond(bond_type="harmonic", + bond_parameters=self.harmonic_bond_params, + particle_pairs=[['A', 'B'], ['B', 'C']]) + pmb.define_angle(angle_type="harmonic", + angle_parameters=self.harmonic_angle_params, + particle_triplets=[('A', 'B', 'C')]) + pmb.define_residue(name="R_A", + central_bead="A", + side_chains=[]) + pmb.define_residue(name="R_B", + central_bead="B", + side_chains=[]) + pmb.define_residue(name="R_C", + central_bead="C", + side_chains=[]) + pmb.define_molecule(name="M_angle", + residue_list=["R_A", "R_B", "R_C"]) + + molecule_ids = pmb.create_molecule(name="M_angle", + number_of_molecules=1, + espresso_system=espresso_system, + gen_angle=True) + + particle_ids = pmb.db._find_instance_ids_by_attribute(pmb_type="particle", + attribute="molecule_id", + value=molecule_ids[0]) + central_particle_id = next(pid for pid in particle_ids + if pmb.db.get_instance(pmb_type="particle", + instance_id=pid).name == "B") + + angle_object = self.get_angle_object(central_particle_id=central_particle_id) + self.assertIsNotNone(angle_object, "No angle object found on the central particle") + self.assertEqual(len(pmb.get_instances_df(pmb_type="angle")), 1) + + pmb.delete_instances_in_system(instance_id=molecule_ids[0], + pmb_type="molecule", + espresso_system=espresso_system) + pmb.db.delete_templates(pmb_type="angle") + pmb.db.delete_templates(pmb_type="bond") + if __name__ == '__main__': ut.main() From 37266756186d6a55e1626aa958c6351847b8a45e Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Tue, 21 Apr 2026 16:15:37 +0200 Subject: [PATCH 04/13] in/output of angle pmb_type --- pyMBE/storage/io.py | 48 ++++++++++++++++++++++++++++++-- pyMBE/storage/manager.py | 2 +- testsuite/database_unit_tests.py | 24 +++++++++++++++- testsuite/test_io_database.py | 26 ++++++++++++++++- 4 files changed, 95 insertions(+), 5 deletions(-) diff --git a/pyMBE/storage/io.py b/pyMBE/storage/io.py index 70826528..ec97d7d6 100644 --- a/pyMBE/storage/io.py +++ b/pyMBE/storage/io.py @@ -29,10 +29,12 @@ from pyMBE.storage.templates.residue import ResidueTemplate from pyMBE.storage.templates.molecule import MoleculeTemplate from pyMBE.storage.templates.bond import BondTemplate +from pyMBE.storage.templates.angle import AngleTemplate from pyMBE.storage.instances.particle import ParticleInstance from pyMBE.storage.instances.residue import ResidueInstance from pyMBE.storage.instances.molecule import MoleculeInstance from pyMBE.storage.instances.bond import BondInstance +from pyMBE.storage.instances.angle import AngleInstance from pyMBE.storage.reactions.reaction import Reaction, ReactionParticipant from pyMBE.storage.templates.peptide import PeptideTemplate from pyMBE.storage.instances.peptide import PeptideInstance @@ -123,7 +125,7 @@ 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, angle, and hydrogel types. """ folder = Path(folder) if not folder.exists(): @@ -134,6 +136,7 @@ def _load_database_csv(db, folder): "residue", "molecule", "bond", + "angle", "peptide", "protein", "hydrogel", @@ -218,6 +221,22 @@ def _load_database_csv(db, folder): particle_name2=None if particle_name2 == "" else particle_name2, parameters=parameters) templates[tpl.name] = tpl + elif pmb_type == "angle": + params_raw = _decode(row.get("parameters", "")) or {} + parameters: Dict[str, Any] = {} + side_particle1 = row.get("side_particle1", "") or "" + central_particle = row.get("central_particle", "") or "" + side_particle2 = row.get("side_particle2", "") or "" + for k, v in params_raw.items(): + if isinstance(v, dict) and {"magnitude", "units", "dimension"}.issubset(v.keys()): + parameters[k] = PintQuantity.from_dict(v) + tpl = AngleTemplate(name=row["name"], + angle_type=row.get("angle_type", ""), + side_particle1=None if side_particle1 == "" else side_particle1, + central_particle=None if central_particle == "" else central_particle, + side_particle2=None if side_particle2 == "" else side_particle2, + parameters=parameters) + templates[tpl.name] = tpl elif pmb_type == "hydrogel": node_map_raw = _decode(row.get("node_map", "")) or [] chain_map_raw = _decode(row.get("chain_map", "")) or [] @@ -308,6 +327,13 @@ def _load_database_csv(db, folder): particle_id1=int(row["particle_id1"]), particle_id2=int(row["particle_id2"])) instances[inst.bond_id] = inst + elif pmb_type == "angle": + inst = AngleInstance(name=row["name"], + angle_id=int(row["angle_id"]), + particle_id1=int(row["particle_id1"]), + particle_id2=int(row["particle_id2"]), + particle_id3=int(row["particle_id3"])) + instances[inst.angle_id] = inst elif pmb_type == "hydrogel": inst = HydrogelInstance(name=row["name"], assembly_id=int(row["assembly_id"])) @@ -401,6 +427,17 @@ def _save_database_csv(db, folder): "particle_name2": tpl.particle_name2, "bond_type": tpl.bond_type, "parameters": _encode(params_serial)}) + elif pmb_type == "angle" and isinstance(tpl, AngleTemplate): + params_serial = {} + for k, v in tpl.parameters.items(): + if isinstance(v, PintQuantity): + params_serial[k] = v.to_dict() + rows.append({"name": tpl.name, + "side_particle1": tpl.side_particle1, + "central_particle": tpl.central_particle, + "side_particle2": tpl.side_particle2, + "angle_type": tpl.angle_type, + "parameters": _encode(params_serial)}) # HYDROGEL TEMPLATE elif pmb_type == "hydrogel" and isinstance(tpl, HydrogelTemplate): rows.append({"name": tpl.name, @@ -465,6 +502,13 @@ def _save_database_csv(db, folder): "bond_id": int(inst.bond_id), "particle_id1": int(inst.particle_id1), "particle_id2": int(inst.particle_id2)}) + elif pmb_type == "angle" and isinstance(inst, AngleInstance): + rows.append({"pmb_type": pmb_type, + "name": inst.name, + "angle_id": int(inst.angle_id), + "particle_id1": int(inst.particle_id1), + "particle_id2": int(inst.particle_id2), + "particle_id3": int(inst.particle_id3)}) elif pmb_type == "hydrogel" and isinstance(inst, HydrogelInstance): rows.append({"pmb_type": pmb_type, "name": inst.name, @@ -488,4 +532,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 8131303e..f8ad826a 100644 --- a/pyMBE/storage/manager.py +++ b/pyMBE/storage/manager.py @@ -90,7 +90,7 @@ def __init__(self,units): "peptide", "protein"] self._assembly_like_types = ["hydrogel"] - self._pmb_types = ["particle", "residue"] + self._molecule_like_types + self._assembly_like_types + self._pmb_types = ["particle", "residue", "angle"] + self._molecule_like_types + self._assembly_like_types self.espresso_bond_instances= {} self.espresso_angle_instances= {} diff --git a/testsuite/database_unit_tests.py b/testsuite/database_unit_tests.py index af3f3b10..46e01a15 100644 --- a/testsuite/database_unit_tests.py +++ b/testsuite/database_unit_tests.py @@ -26,8 +26,10 @@ from pyMBE.storage.instances.peptide import PeptideInstance from pyMBE.storage.instances.protein import ProteinInstance from pyMBE.storage.instances.bond import BondInstance +from pyMBE.storage.instances.angle import AngleInstance from pyMBE.storage.instances.hydrogel import HydrogelInstance from pyMBE.storage.templates.bond import BondTemplate +from pyMBE.storage.templates.angle import AngleTemplate from pyMBE.storage.templates.hydrogel import HydrogelNode from pyMBE.storage.pint_quantity import PintQuantity from pyMBE.storage.reactions.reaction import Reaction, ReactionParticipant @@ -410,6 +412,14 @@ def test_instance_id_validators(self): self.assertRaises(ValueError, BondInstance, **inputs) + inputs = {"name":"A", + "angle_id":-1, + "particle_id1":1, + "particle_id2":2, + "particle_id3":3} + self.assertRaises(ValueError, + AngleInstance, + **inputs) def test_make_name_bond_template(self): inputs = {"bond_type": "harmonic", @@ -419,6 +429,18 @@ def test_make_name_bond_template(self): bond_tpl = BondTemplate(**inputs) self.assertRaises(RuntimeError, bond_tpl._make_name) + + def test_make_name_angle_template(self): + inputs = {"angle_type": "harmonic", + "parameters": {"k": PintQuantity(magnitude=1, + units="reduced_energy", + dimension="energy"), + "phi_0": PintQuantity(magnitude=1, + units="dimensionless", + dimension="dimensionless")}} + angle_tpl = AngleTemplate(**inputs) + self.assertRaises(RuntimeError, + angle_tpl._make_name) def test_exceptions_pint_quantity(self): units = pint.UnitRegistry() @@ -480,4 +502,4 @@ def test_exceptions_reaction_template(self): **inputs) if __name__ == '__main__': - ut.main() \ No newline at end of file + ut.main() diff --git a/testsuite/test_io_database.py b/testsuite/test_io_database.py index 7c03c522..80c3ba34 100644 --- a/testsuite/test_io_database.py +++ b/testsuite/test_io_database.py @@ -293,6 +293,31 @@ def test_io_bond_templates(self): pd.testing.assert_frame_equal(pmb.get_templates_df(pmb_type="bond"), new_pmb.get_templates_df(pmb_type="bond")) + def test_io_angle_templates(self): + """ + Checks that information in the pyMBE database about + angle templates is stored to file and loaded back properly. + """ + pmb = pyMBE.pymbe_library(seed=42) + parameters1 = {"k": 50.0 * pmb.units.reduced_energy, + "phi_0": 1.0 * pmb.units("")} + parameters2 = {"k": 30.0 * pmb.units.reduced_energy, + "phi_0": 0.5 * pmb.units("")} + pmb.define_angle(angle_type="harmonic", + angle_parameters=parameters1, + particle_triplets=[("Z", "X", "Z"), + ("X", "Z", "Y")]) + pmb.define_default_angle(angle_type="cosine", + angle_parameters=parameters2) + new_pmb = pyMBE.pymbe_library(23) + with tempfile.TemporaryDirectory() as tmp_directory: + # Save and load the database + pmb.save_database(tmp_directory) + new_pmb.load_database(tmp_directory) + # Test that the loaded angle templates are equal to the original + pd.testing.assert_frame_equal(pmb.get_templates_df(pmb_type="angle"), + new_pmb.get_templates_df(pmb_type="angle")) + def test_io_residues_templates(self): """ Checks that information in the pyMBE database about @@ -755,4 +780,3 @@ def test_load_database_missing_folder(self): if __name__ == '__main__': ut.main() - From e2f950dea8a26ba73436c0a9fd251c39008b284e Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Wed, 22 Apr 2026 13:42:41 +0200 Subject: [PATCH 05/13] uncovered io.py fixed --- testsuite/test_io_database.py | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/testsuite/test_io_database.py b/testsuite/test_io_database.py index 80c3ba34..fcbd7b84 100644 --- a/testsuite/test_io_database.py +++ b/testsuite/test_io_database.py @@ -462,7 +462,8 @@ def test_io_instances(self): instances created into espresso is stored to file and loaded back properly. """ pmb = pyMBE.pymbe_library(seed=42) - # Test instances of a hydrogel (tests hydrogel, molecule, residue, bond and particle instances) + # Test instances of a hydrogel plus a standalone residue with an angle + # (tests hydrogel, molecule, residue, bond, angle and particle instances) pmb.define_particle(name="Z", sigma=3.5 * pmb.units.reduced_length, cutoff=4 * pmb.units.reduced_length, @@ -486,8 +487,16 @@ def test_io_instances(self): particle_pairs=[["Z","Z"], ["Z","X"], ["X","X"]]) + angle_parameters = {"k": 50.0 * pmb.units.reduced_energy, + "phi_0": 1.0 * pmb.units("")} + pmb.define_angle(angle_type="harmonic", + angle_parameters=angle_parameters, + particle_triplets=[("X", "Z", "Z")]) pmb.define_molecule(name="M1", residue_list=["R1"]*1) + angle_residue_id = pmb.create_residue(name="R1", + espresso_system=espresso_system, + gen_angle=True) diamond_lattice = DiamondLattice(4, 3.5 * pmb.units.reduced_length) lattice_builder = pmb.initialize_lattice_builder(diamond_lattice) # Setting up node topology --> Nodes are particles of type "X" @@ -522,16 +531,22 @@ def test_io_instances(self): new_pmb.get_instances_df(pmb_type="hydrogel")) pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="molecule"), new_pmb.get_instances_df(pmb_type="molecule")) - pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="residues"), - new_pmb.get_instances_df(pmb_type="residues")) + pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="residue"), + new_pmb.get_instances_df(pmb_type="residue")) pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="bond"), new_pmb.get_instances_df(pmb_type="bond")) + pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="angle"), + new_pmb.get_instances_df(pmb_type="angle")) pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="particle"), new_pmb.get_instances_df(pmb_type="particle")) # Clean up before the next test pmb.delete_instances_in_system(espresso_system=espresso_system, instance_id=assembly_id, pmb_type="hydrogel") + pmb.delete_instances_in_system(espresso_system=espresso_system, + instance_id=angle_residue_id, + pmb_type="residue") + pmb.db.delete_templates(pmb_type="angle") # Test instances of a peptide (tests peptide, residue, bond and particle instances) path_to_interactions=pmb.root / "parameters" / "peptides" / "Lunkad2021" path_to_pka=pmb.root / "parameters" / "pka_sets" / "Hass2015.json" @@ -564,8 +579,8 @@ def test_io_instances(self): # Test that the loaded instances are equal to the original pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="peptide"), new_pmb.get_instances_df(pmb_type="peptide")) - pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="residues"), - new_pmb.get_instances_df(pmb_type="residues")) + pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="residue"), + new_pmb.get_instances_df(pmb_type="residue")) pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="bond"), new_pmb.get_instances_df(pmb_type="bond")) pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="particle"), @@ -605,8 +620,8 @@ def test_io_instances(self): # Test that the loaded instances are equal to the original pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="protein"), new_pmb.get_instances_df(pmb_type="protein")) - pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="residues"), - new_pmb.get_instances_df(pmb_type="residues")) + pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="residue"), + new_pmb.get_instances_df(pmb_type="residue")) pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="bond"), new_pmb.get_instances_df(pmb_type="bond")) pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="particle"), From 28f1e4a0ec74ef799a6c8fe2c68dcfa4f30789ec Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sat, 2 May 2026 04:07:46 +0200 Subject: [PATCH 06/13] angle potential support for hydrogels --- pyMBE/pyMBE.py | 149 ++++++++++++++--- samples/build_hydrogel_angle.py | 161 +++++++++++++++++++ testsuite/CTestTestfile.cmake | 1 + testsuite/hydrogel_builder.py | 2 - testsuite/hydrogel_builder_with_angles.py | 186 ++++++++++++++++++++++ testsuite/lattice_builder.py | 47 +++++- 6 files changed, 520 insertions(+), 26 deletions(-) create mode 100644 samples/build_hydrogel_angle.py create mode 100644 testsuite/hydrogel_builder_with_angles.py diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index f91dd244..0636eb31 100644 --- a/pyMBE/pyMBE.py +++ b/pyMBE/pyMBE.py @@ -227,7 +227,7 @@ def _create_espresso_bond_instance(self, bond_type, bond_parameters): d_r_max = bond_parameters["d_r_max"].m_as("reduced_length")) return bond_instance - def _create_hydrogel_chain(self, hydrogel_chain, nodes, espresso_system, use_default_bond=False): + def _create_hydrogel_chain(self, hydrogel_chain, nodes, espresso_system, use_default_bond=False, gen_angle=False): """ Creates a chain between two nodes of a hydrogel. @@ -243,6 +243,11 @@ def _create_hydrogel_chain(self, hydrogel_chain, nodes, espresso_system, use_def use_default_bond ('bool', optional): If True, use a default bond template if no specific template exists. Defaults to False. + gen_angle ('bool', optional): + If True, generate the angle potentials internal to the created + chain molecule. Junction angles near the hydrogel crosslinkers + are handled separately at the hydrogel level. + Return: ('int'): molecule_id of the created hydrogel chian. @@ -262,8 +267,8 @@ def _create_hydrogel_chain(self, hydrogel_chain, nodes, espresso_system, use_def node_start_label = self.lattice_builder._create_node_label(node_start) node_end_label = self.lattice_builder._create_node_label(node_end) _, reverse = self.lattice_builder._get_node_vector_pair(node_start, node_end) - if node_start != node_end or residue_list == residue_list[::-1]: - ValueError(f"Aborted creation of hydrogel chain between '{node_start}' and '{node_end}' because pyMBE could not resolve a unique topology for that chain") + if node_start == node_end and residue_list != residue_list[::-1]: + raise ValueError(f"Aborted creation of hydrogel chain between '{node_start}' and '{node_end}' because pyMBE could not resolve a unique topology for that chain") if reverse: reverse_residue_order=True else: @@ -298,25 +303,83 @@ def _create_hydrogel_chain(self, hydrogel_chain, nodes, espresso_system, use_def list_of_first_residue_positions=[first_bead_pos.tolist()], #Start at the first node backbone_vector=np.array(backbone_vector)/l0, use_default_bond=use_default_bond, - reverse_residue_order=reverse_residue_order)[0] + reverse_residue_order=reverse_residue_order, + gen_angle=gen_angle)[0] # Bond chain to the hydrogel nodes chain_pids = self.db._find_instance_ids_by_attribute(pmb_type="particle", attribute="molecule_id", value=mol_id) - bond_tpl1 = self.get_bond_template(particle_name1=nodes[node_start_label]["name"], - particle_name2=part_start_chain_name, - use_default_bond=use_default_bond) - start_bond_instance = self._get_espresso_bond_instance(bond_template=bond_tpl1, - espresso_system=espresso_system) - bond_tpl2 = self.get_bond_template(particle_name1=nodes[node_end_label]["name"], - particle_name2=part_end_chain_name, - use_default_bond=use_default_bond) - end_bond_instance = self._get_espresso_bond_instance(bond_template=bond_tpl2, - espresso_system=espresso_system) - espresso_system.part.by_id(start_node_id).add_bond((start_bond_instance, chain_pids[0])) - espresso_system.part.by_id(chain_pids[-1]).add_bond((end_bond_instance, end_node_id)) + self.create_bond(particle_id1=start_node_id, + particle_id2=chain_pids[0], + espresso_system=espresso_system, + use_default_bond=use_default_bond) + self.create_bond(particle_id1=chain_pids[-1], + particle_id2=end_node_id, + espresso_system=espresso_system, + use_default_bond=use_default_bond) return mol_id + def _generate_hydrogel_crosslinker_angles(self, espresso_system, central_particle_ids): + """ + Generate hydrogel angles centered on crosslinkers and adjacent terminal beads. + + If the user defines any explicit angle template for such junction + triplets, then all required junction triplets must be defined. If none + are defined, hydrogel construction proceeds without crosslinker-adjacent + angles. + """ + particle_instances = self.db.get_instances(pmb_type="particle") + bonded_neighbors = {} + for bond in self.db.get_instances(pmb_type="bond").values(): + bonded_neighbors.setdefault(bond.particle_id1, set()).add(bond.particle_id2) + bonded_neighbors.setdefault(bond.particle_id2, set()).add(bond.particle_id1) + + triplets = [] + for central_particle_id in sorted(set(central_particle_ids)): + neighbors = sorted(bonded_neighbors.get(central_particle_id, set())) + central_name = particle_instances[central_particle_id].name + for idx_i in range(len(neighbors)): + for idx_k in range(idx_i + 1, len(neighbors)): + side_particle_id1 = neighbors[idx_i] + side_particle_id3 = neighbors[idx_k] + side_name1 = particle_instances[side_particle_id1].name + side_name3 = particle_instances[side_particle_id3].name + angle_key = AngleTemplate.make_angle_key(side1=side_name1, + central=central_name, + side2=side_name3) + triplets.append((side_particle_id1, + central_particle_id, + side_particle_id3, + angle_key)) + + defined_angle_templates = self.db.get_templates(pmb_type="angle") + defined_angle_keys = { + angle_key + for _, _, _, angle_key in triplets + if angle_key in defined_angle_templates + } + if not defined_angle_keys: + logging.warning("No angle templates defined for hydrogel crosslinkers") + return + + missing_angle_keys = sorted({ + angle_key + for _, _, _, angle_key in triplets + if angle_key not in defined_angle_keys + }) + if missing_angle_keys: + raise ValueError( + "Hydrogel crosslinker-adjacent angle templates must be defined for all required triplets. " + f"Missing definitions for: {missing_angle_keys}" + ) + + for side_particle_id1, central_particle_id, side_particle_id3, _ in triplets: + self.create_angle(particle_id1=side_particle_id1, + particle_id2=central_particle_id, + particle_id3=side_particle_id3, + espresso_system=espresso_system, + use_default_angle=False) + def _create_hydrogel_node(self, node_index, node_name, espresso_system): """ Set a node residue type. @@ -903,7 +966,7 @@ def create_counterions(self, object_name, cation_name, anion_name, espresso_syst logging.info(f'Ion type: {name} created number: {counterion_number[name]}') return counterion_number - def create_hydrogel(self, name, espresso_system, use_default_bond=False): + def create_hydrogel(self, name, espresso_system, use_default_bond=False, gen_angle=False): """ Creates a hydrogel in espresso_system using a pyMBE hydrogel template given by 'name' @@ -917,6 +980,11 @@ def create_hydrogel(self, name, espresso_system, use_default_bond=False): use_default_bond ('bool', optional): If True, use a default bond template if no specific template exists. Defaults to False. + gen_angle ('bool', optional): + If True, generate angle potentials for the internal hydrogel + chains and, when explicitly defined, for all crosslinker-adjacent + triplets. Defaults to False. + Returns: ('int'): id of the hydrogel instance created. """ @@ -927,6 +995,7 @@ def create_hydrogel(self, name, espresso_system, use_default_bond=False): assembly_id = self.db._propose_instance_id(pmb_type="hydrogel") # Create the nodes nodes = {} + hydrogel_angle_centers = set() node_topology = hydrogel_tpl.node_map for node in node_topology: node_index = node.lattice_index @@ -944,15 +1013,59 @@ def create_hydrogel(self, name, espresso_system, use_default_bond=False): molecule_id = self._create_hydrogel_chain(hydrogel_chain=hydrogel_chain, nodes=nodes, espresso_system=espresso_system, - use_default_bond=use_default_bond) + use_default_bond=use_default_bond, + gen_angle=gen_angle) self.db._update_instance(instance_id=molecule_id, pmb_type="molecule", attribute="assembly_id", value=assembly_id) + if gen_angle: + residue_ids = self.db._find_instance_ids_by_attribute(pmb_type="residue", + attribute="molecule_id", + value=molecule_id) + first_residue_id = min(residue_ids) + last_residue_id = max(residue_ids) + first_residue = self.db.get_instance(pmb_type="residue", + instance_id=first_residue_id) + last_residue = self.db.get_instance(pmb_type="residue", + instance_id=last_residue_id) + first_central_bead_name = self.db.get_template(pmb_type="residue", + name=first_residue.name).central_bead + last_central_bead_name = self.db.get_template(pmb_type="residue", + name=last_residue.name).central_bead + particle_instances = self.db.get_instances(pmb_type="particle") + first_residue_particle_ids = self.db._find_instance_ids_by_attribute(pmb_type="particle", + attribute="residue_id", + value=first_residue_id) + last_residue_particle_ids = self.db._find_instance_ids_by_attribute(pmb_type="particle", + attribute="residue_id", + value=last_residue_id) + first_bead_id = None + for particle_id in first_residue_particle_ids: + if particle_instances[particle_id].name == first_central_bead_name: + first_bead_id = particle_id + break + + last_bead_id = None + for particle_id in last_residue_particle_ids: + if particle_instances[particle_id].name == last_central_bead_name: + last_bead_id = particle_id + break + node_start_label = self.lattice_builder._create_node_label(hydrogel_chain.node_start) + node_end_label = self.lattice_builder._create_node_label(hydrogel_chain.node_end) + hydrogel_angle_centers.update({ + nodes[node_start_label]["id"], + nodes[node_end_label]["id"], + first_bead_id, + last_bead_id, + }) self.db._propagate_id(root_type="hydrogel", root_id=assembly_id, attribute="assembly_id", value=assembly_id) + if gen_angle: + self._generate_hydrogel_crosslinker_angles(espresso_system=espresso_system, + central_particle_ids=hydrogel_angle_centers) # Register an hydrogel instance in the pyMBE databasegit self.db._register_instance(HydrogelInstance(name=name, assembly_id=assembly_id)) diff --git a/samples/build_hydrogel_angle.py b/samples/build_hydrogel_angle.py new file mode 100644 index 00000000..bae855ec --- /dev/null +++ b/samples/build_hydrogel_angle.py @@ -0,0 +1,161 @@ +# +# Example: demonstrating hydrogel angle generation in pyMBE. +# +# This script builds a simple hydrogel of: +# - N: crosslinker/node particles +# - C: residue/chain particles +# +# By default, it defines only C-C-C and generates only chain-internal angles. +# If run with `--include`, it also defines the crosslinker-adjacent +# angle templates N-C-C and C-N-C. +# +# Note: +# - The canonical angle name for N-C-C is C-C-N, because pyMBE sorts the +# two side particles alphabetically when building angle template names. +# + +import argparse +import numpy as np +import espressomd + +import pyMBE +from pyMBE.lib.lattice import DiamondLattice + + +def define_simple_hydrogel_system(include_crosslinker_angles): + """ + Build a simple hydrogel with C residues and N crosslinkers. + """ + pmb = pyMBE.pymbe_library(seed=42) + + node_type = "N" + bead_type = "C" + residue_name = "Res" + molecule_name = "hydrogel_chain" + + bond_length = 0.355 * pmb.units.nm + mpc = 4 + + pmb.define_particle(name=node_type, + sigma=bond_length, + epsilon=1 * pmb.units("reduced_energy")) + pmb.define_particle(name=bead_type, + sigma=bond_length, + epsilon=1 * pmb.units("reduced_energy")) + + pmb.define_residue(name=residue_name, + central_bead=bead_type, + side_chains=[]) + pmb.define_molecule(name=molecule_name, + residue_list=[residue_name] * mpc) + + harmonic_bond_parameters = { + "r_0": bond_length, + "k": 400 * pmb.units("reduced_energy / reduced_length**2"), + } + pmb.define_bond(bond_type="harmonic", + bond_parameters=harmonic_bond_parameters, + particle_pairs=[ + [bead_type, bead_type], + [node_type, bead_type], + ]) + + angle_parameters = { + "k": 25 * pmb.units("reduced_energy"), + "phi_0": np.pi * pmb.units(""), + } + pmb.define_angle(angle_type="harmonic", + angle_parameters=angle_parameters, + particle_triplets=[(bead_type, bead_type, bead_type)]) + + if include_crosslinker_angles: + pmb.define_angle(angle_type="harmonic", + angle_parameters=angle_parameters, + particle_triplets=[ + (node_type, bead_type, bead_type), # canonical: C-C-N + (bead_type, node_type, bead_type), # canonical: C-N-C + ]) + + diamond_lattice = DiamondLattice(mpc, bond_length) + espresso_system = espressomd.System(box_l=[diamond_lattice.box_l] * 3) + pmb.initialize_lattice_builder(diamond_lattice) + + node_topology = [ + {"particle_name": node_type, "lattice_index": index} + for index in diamond_lattice.indices + ] + chain_topology = [] + for node_connectivity in diamond_lattice.connectivity: + node_start = str(diamond_lattice.indices[node_connectivity[0]]) + node_end = str(diamond_lattice.indices[node_connectivity[1]]) + chain_topology.append({ + "node_start": node_start, + "node_end": node_end, + "molecule_name": molecule_name, + }) + + pmb.define_hydrogel(name="simple_hydrogel", + node_map=node_topology, + chain_map=chain_topology) + hydrogel_id = pmb.create_hydrogel(name="simple_hydrogel", + espresso_system=espresso_system, + gen_angle=True) + return pmb, espresso_system, hydrogel_id + + +def summarize_angles(pmb): + """ + Return angle templates, angle instances, and counts by canonical angle name. + """ + angle_templates_df = pmb.get_templates_df(pmb_type="angle") + angle_instances_df = pmb.get_instances_df(pmb_type="angle") + angle_counts = angle_instances_df["name"].value_counts().sort_index() + return angle_templates_df, angle_instances_df, angle_counts + + +def run_case(label, include_crosslinker_angles, expected_angle_names): + """ + Build one hydrogel case, print angle counts, and assert expectations. + """ + print(f"\n############ {label} ############") + pmb, _, hydrogel_id = define_simple_hydrogel_system( + include_crosslinker_angles=include_crosslinker_angles + ) + angle_templates_df, angle_instances_df, angle_counts = summarize_angles(pmb) + print(f"Hydrogel assembly id: {hydrogel_id}") + print("\nAngle templates:") + print(angle_templates_df) + print("\nAngle instance dataframe:") + print(angle_instances_df) + print("Angle instance counts:") + print(angle_counts) + + angle_names = set(angle_counts.index) + if angle_names != expected_angle_names: + raise AssertionError( + f"{label}: expected angle names {sorted(expected_angle_names)}, " + f"got {sorted(angle_names)}" + ) + + print(f"{label}: OK") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Build a simple hydrogel and optionally include crosslinker-adjacent angles." + ) + parser.add_argument( + "--include", + action="store_true", + help="Include crosslinker-adjacent angles (N-C-C and C-N-C).", + ) + args = parser.parse_args() + + if args.include: + run_case(label="Chain + crosslinker angles", + include_crosslinker_angles=True, + expected_angle_names={"C-C-C", "C-C-N", "C-N-C"}) + else: + run_case(label="Chain angles only", + include_crosslinker_angles=False, + expected_angle_names={"C-C-C"}) diff --git a/testsuite/CTestTestfile.cmake b/testsuite/CTestTestfile.cmake index 4dfeea96..63caeaea 100644 --- a/testsuite/CTestTestfile.cmake +++ b/testsuite/CTestTestfile.cmake @@ -76,4 +76,5 @@ pymbe_add_test(PATH determine_reservoir_concentrations_unit_test.py) 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 hydrogel_builder_with_angles.py) pymbe_add_test(PATH database_unit_tests.py) diff --git a/testsuite/hydrogel_builder.py b/testsuite/hydrogel_builder.py index eeb0fab5..da4017cc 100644 --- a/testsuite/hydrogel_builder.py +++ b/testsuite/hydrogel_builder.py @@ -111,8 +111,6 @@ name=hydrogel_name) hydrogel_inst = pmb.db.get_instance(pmb_type="hydrogel", instance_id=hydrogel_id) - - class Test(ut.TestCase): def test_create_hydrogel_missing_template(self): """ diff --git a/testsuite/hydrogel_builder_with_angles.py b/testsuite/hydrogel_builder_with_angles.py new file mode 100644 index 00000000..658c89d5 --- /dev/null +++ b/testsuite/hydrogel_builder_with_angles.py @@ -0,0 +1,186 @@ +# +# Copyright (C) 2024-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 collections import Counter +import numpy as np +import unittest as ut +import pyMBE +from pyMBE.lib.lattice import DiamondLattice +import espressomd + +espresso_system = espressomd.System(box_l=[10] * 3) + + +def build_simple_hydrogel_with_optional_angles(junction_angle_mode="none", mpc_local=4): + """ + Build a simple hydrogel used to validate hydrogel-specific angle support. + """ + pmb_local = pyMBE.pymbe_library(seed=7) + node_type = "N" + bead_type = "C" + residue_name = "Res" + molecule_name = "simple_chain" + bond_length = 0.355 * pmb_local.units.nm + harmonic_bond_parameters = { + "r_0": bond_length, + "k": 400 * pmb_local.units("reduced_energy / reduced_length**2"), + } + + pmb_local.define_particle(name=node_type, + sigma=bond_length, + epsilon=1 * pmb_local.units("reduced_energy")) + pmb_local.define_particle(name=bead_type, + sigma=bond_length, + epsilon=1 * pmb_local.units("reduced_energy")) + pmb_local.define_residue(name=residue_name, + central_bead=bead_type, + side_chains=[]) + pmb_local.define_molecule(name=molecule_name, + residue_list=[residue_name] * mpc_local) + pmb_local.define_bond(bond_type="harmonic", + bond_parameters=harmonic_bond_parameters, + particle_pairs=[[bead_type, bead_type], + [bead_type, node_type]]) + + angle_parameters = { + "k": 5 * pmb_local.units("reduced_energy"), + "phi_0": np.pi * pmb_local.units(""), + } + pmb_local.define_angle(angle_type="harmonic", + angle_parameters=angle_parameters, + particle_triplets=[(bead_type, bead_type, bead_type)]) + + if junction_angle_mode in {"partial", "full"}: + particle_triplets = [(bead_type, node_type, bead_type)] + if junction_angle_mode == "full": + particle_triplets.append((node_type, bead_type, bead_type)) + pmb_local.define_angle(angle_type="harmonic", + angle_parameters=angle_parameters, + particle_triplets=particle_triplets) + + diamond_lattice_local = DiamondLattice(mpc_local, bond_length) + espresso_system.box_l = [diamond_lattice_local.box_l] * 3 + pmb_local.initialize_lattice_builder(diamond_lattice_local) + + node_topology_local = [ + {"particle_name": node_type, "lattice_index": index} + for index in diamond_lattice_local.indices + ] + chain_topology_local = [] + for node_connectivity in diamond_lattice_local.connectivity: + node_start = str(diamond_lattice_local.indices[node_connectivity[0]]) + node_end = str(diamond_lattice_local.indices[node_connectivity[1]]) + chain_topology_local.append({ + "node_start": node_start, + "node_end": node_end, + "molecule_name": molecule_name, + }) + + hydrogel_name_local = "simple_hydrogel" + pmb_local.define_hydrogel(hydrogel_name_local, + node_topology_local, + chain_topology_local) + return pmb_local, espresso_system, hydrogel_name_local, diamond_lattice_local + + +def get_angle_counts(pmb_local): + """ + Return counts of angle instances by canonical angle name. + """ + return Counter( + angle.name + for angle in pmb_local.db.get_instances(pmb_type="angle").values() + ) + + +def expected_crosslinker_angle_counts(diamond_lattice_local): + """ + Return expected chain and crosslinker angle counts for the simple hydrogel. + """ + number_of_nodes = 8 + number_of_chains = 16 + crosslinker_functionality = 4 + + return { + "C-C-C": number_of_chains * (diamond_lattice_local.mpc - 2), + "C-C-N": 2 * number_of_chains, + "C-N-C": number_of_nodes * (crosslinker_functionality * (crosslinker_functionality - 1) // 2), + } + + +class Test(ut.TestCase): + def setUp(self): + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + espresso_system.non_bonded_inter.reset() + + def test_hydrogel_gen_angle_defaults_to_false(self): + """ + Hydrogel creation should not generate angles unless gen_angle is requested. + """ + pmb_local, espresso_system_local, hydrogel_name_local, _ = build_simple_hydrogel_with_optional_angles( + junction_angle_mode="full" + ) + pmb_local.create_hydrogel(hydrogel_name_local, espresso_system_local) + self.assertEqual(len(pmb_local.db.get_instances(pmb_type="angle")), 0) + + def test_hydrogel_chain_angles_are_created_without_crosslinker_templates(self): + """ + If only chain angles are defined, hydrogel creation should still succeed + and create only the chain-internal angles. + """ + pmb_local, espresso_system_local, hydrogel_name_local, diamond_lattice_local = build_simple_hydrogel_with_optional_angles( + junction_angle_mode="none" + ) + with self.assertLogs(level="WARNING") as log_context: + pmb_local.create_hydrogel(hydrogel_name_local, espresso_system_local, gen_angle=True) + + angle_counts = get_angle_counts(pmb_local) + self.assertEqual( + angle_counts, + {"C-C-C": len(diamond_lattice_local.connectivity) * (diamond_lattice_local.mpc - 2)}, + ) + self.assertTrue( + any("No angle templates defined for hydrogel crosslinkers" in message for message in log_context.output) + ) + + def test_hydrogel_junction_angle_counts_are_correct_when_defined(self): + """ + Hydrogel creation should add the correct number of junction angles when + the explicit templates exist. + """ + pmb_local, espresso_system_local, hydrogel_name_local, diamond_lattice_local = build_simple_hydrogel_with_optional_angles( + junction_angle_mode="full" + ) + pmb_local.create_hydrogel(hydrogel_name_local, espresso_system_local, gen_angle=True) + self.assertEqual(get_angle_counts(pmb_local), + expected_crosslinker_angle_counts(diamond_lattice_local)) + + def test_hydrogel_partial_crosslinker_angle_definitions_raise(self): + """ + Defining only a subset of required crosslinker-adjacent angles should fail. + """ + pmb_local, espresso_system_local, hydrogel_name_local, _ = build_simple_hydrogel_with_optional_angles( + junction_angle_mode="partial" + ) + with self.assertRaises(ValueError): + pmb_local.create_hydrogel(hydrogel_name_local, espresso_system_local, gen_angle=True) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/lattice_builder.py b/testsuite/lattice_builder.py index b84fdeb0..2788ab7a 100644 --- a/testsuite/lattice_builder.py +++ b/testsuite/lattice_builder.py @@ -226,13 +226,17 @@ def test_lattice_setup(self): "k": 400 * pmb2.units('reduced_energy / reduced_length**2')}) # Define nodes dictionary as expected by _create_hydrogel_chain + # using actual hydrogel node particles registered in pyMBE. nodes = {} - id = 0 for label, index in lattice.node_labels.items(): - nodes[label] = {"name": lattice.get_node(label), - "pos": lattice.lattice.indices[index], - "id": id} - id +=1 + node_index = lattice.lattice.indices[index] + node_name = lattice.get_node(label) + node_pos, node_id = pmb2._create_hydrogel_node(node_index=node_index, + node_name=node_name, + espresso_system=espresso_system) + nodes[label] = {"name": node_name, + "pos": node_pos, + "id": node_id} from pyMBE.storage.templates.hydrogel import HydrogelChain # Define hydrogel chain template (reverse geometry) hydrogel_chain = HydrogelChain(node_start=node_b, # reversed on purpose @@ -252,6 +256,38 @@ def test_lattice_setup(self): # Reverse branch MUST reverse the residue list np.testing.assert_equal(actual=created_residues, desired=sequence[::-1], verbose=True) + def test_non_palindromic_self_loop_chain_raises(self): + """ + A self-loop hydrogel chain with a non-palindromic residue sequence is + ambiguous and must be rejected by _create_hydrogel_chain. + """ + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + espresso_system.non_bonded_inter.reset() + + pmb = pyMBE.pymbe_library(seed=11) + define_templates(pmb=pmb) + lattice = pmb.initialize_lattice_builder(diamond) + lattice.strict = False + + non_palindromic_sequence = [Res1, Res2, Res3, Res1] + pmb.define_molecule(name="self_loop_chain", + residue_list=non_palindromic_sequence) + + from pyMBE.storage.templates.hydrogel import HydrogelChain + hydrogel_chain = HydrogelChain(node_start="[0 0 0]", + node_end="[0 0 0]", + molecule_name="self_loop_chain") + + with self.assertRaisesRegex( + ValueError, + "could not resolve a unique topology for that chain", + ): + pmb._create_hydrogel_chain(hydrogel_chain=hydrogel_chain, + nodes={}, + espresso_system=espresso_system, + use_default_bond=False) + def test_plot(self): pmb = pyMBE.pymbe_library(seed=42) diamond = pyMBE.lib.lattice.DiamondLattice(mpc, bond_l) @@ -304,4 +340,3 @@ def test_plot(self): if __name__ == "__main__": ut.main() - From 0732a3ed2890380b2c662dd776a2a4271913691a Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sat, 2 May 2026 04:19:09 +0200 Subject: [PATCH 07/13] removed unused variable --- pyMBE/pyMBE.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index 0636eb31..99a7098e 100644 --- a/pyMBE/pyMBE.py +++ b/pyMBE/pyMBE.py @@ -286,8 +286,6 @@ def _create_hydrogel_chain(self, hydrogel_chain, nodes, espresso_system, use_def name=molecule_name).residue_list part_start_chain_name = self.db.get_template(pmb_type="residue", name=chain_residues[0]).central_bead - part_end_chain_name = self.db.get_template(pmb_type="residue", - name=chain_residues[-1]).central_bead lj_parameters = self.get_lj_parameters(particle_name1=nodes[node_start_label]["name"], particle_name2=part_start_chain_name) bond_tpl = self.get_bond_template(particle_name1=nodes[node_start_label]["name"], From 068e4f27673880a7b7bf23485ce5ed2aa7fe3bb4 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sat, 2 May 2026 19:45:22 +0200 Subject: [PATCH 08/13] copyright corrected --- samples/build_angle.py | 104 ++++++- samples/build_hydrogel.py | 319 ++++++++++++++++------ samples/build_hydrogel_angle.py | 161 ----------- testsuite/angle_tests.py | 2 +- testsuite/hydrogel_builder_with_angles.py | 2 +- 5 files changed, 327 insertions(+), 261 deletions(-) delete mode 100644 samples/build_hydrogel_angle.py diff --git a/samples/build_angle.py b/samples/build_angle.py index 462eff1e..c1465111 100644 --- a/samples/build_angle.py +++ b/samples/build_angle.py @@ -1,17 +1,32 @@ # +# 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 . +# # Example: demonstrating angle potential support in pyMBE. # -# This script shows two complementary usages of `gen_angle=True`: +# This script shows three usages of `gen_angle=True`: # 1. create_residue with gen_angle=True # - builds a single residue with auto-generated angles # 2. create_molecule with gen_angle=True # - builds a chain of residues; angles are auto-generated across the # full molecule, including angles that span residue boundaries -# -# Note: angle auto-generation is currently supported only in `create_residue` -# and `create_molecule`. Other higher-level builders (e.g. `create_hydrogel`) -# do not yet propagate `gen_angle`. -# +# 3. create_residue with a nested residue as side chain +# - builds a residue whose side chain is itself a residue; verifies +# that all bonded triplets (intra- and cross-level) are generated import pyMBE import numpy as np @@ -41,6 +56,11 @@ sigma=0.4 * pmb.units.nm, epsilon=1 * pmb.units('reduced_energy')) +pmb.define_particle(name='D', + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + # ---------------------------------------------------------------------- # Default bond used everywhere # ---------------------------------------------------------------------- @@ -108,11 +128,71 @@ def show_database(label): pmb.define_molecule(name="Mol_1", residue_list=["Res_2"] * 4) -pmb.create_molecule(name="Mol_1", - number_of_molecules=1, - espresso_system=espresso_system, - backbone_vector=[1.0, 0.0, 0.0], - use_default_bond=True, - gen_angle=True) +molecule_ids = pmb.create_molecule(name="Mol_1", + number_of_molecules=1, + espresso_system=espresso_system, + backbone_vector=[1.0, 0.0, 0.0], + use_default_bond=True, + gen_angle=True) show_database("Demo 2: linear molecule with gen_angle=True") + +pmb.delete_instances_in_system(instance_id=molecule_ids[0], + pmb_type="molecule", + espresso_system=espresso_system) + +# ====================================================================== +# Demo 3: create_residue with a nested residue as side chain +# ====================================================================== +# Topology: D₁ where SubRes_3 = B bonded to C and D₂ +# \ +# A +# | +# B -- C +# | +# D₂ +# +# SubRes_3 (central B, side chains C and D) is passed as a side chain of +# the outer residue NestedRes_3 (central A, direct side chain D). +# D₁ and D₂ are distinct instances of particle type D. +# +# All bonded triplets and their canonical angle names: +# centered at A: D₁ -- A -- B → "B-A-D" (outer {B,D} sorted) +# centered at B: A -- B -- C → "A-B-C" +# centered at B: A -- B -- D₂ → "A-B-D" +# centered at B: C -- B -- D₂ → "C-B-D" +pmb.define_residue(name="SubRes_3", + central_bead="B", + side_chains=["C", "D"]) + +pmb.define_residue(name="NestedRes_3", + central_bead="A", + side_chains=["D", "SubRes_3"]) + +nested_residue_id = pmb.create_residue(name="NestedRes_3", + espresso_system=espresso_system, + use_default_bond=True, + gen_angle=True) + +show_database("Demo 3: nested residue with gen_angle=True") + +angle_instances_df = pmb.get_instances_df(pmb_type="angle") +particle_instances_df = pmb.get_instances_df(pmb_type="particle") +print(f"\nAngle instance count: {len(angle_instances_df)}") + +id_to_type = dict(zip(particle_instances_df["particle_id"], + particle_instances_df["name"])) +actual_triplets = set() +for _, row in angle_instances_df.iterrows(): + p1 = id_to_type[row["particle_id1"]] + center = id_to_type[row["particle_id2"]] + p3 = id_to_type[row["particle_id3"]] + outer = sorted([p1, p3]) + actual_triplets.add(f"{outer[0]}-{center}-{outer[1]}") + +print("Resolved angle triplets:", sorted(actual_triplets)) +expected_triplets = {"A-B-C", "A-B-D", "B-A-D", "C-B-D"} +assert actual_triplets == expected_triplets, ( + f"Demo 3: expected {sorted(expected_triplets)}, got {sorted(actual_triplets)}" +) +print("Demo 3: all expected angle triplets generated OK") diff --git a/samples/build_hydrogel.py b/samples/build_hydrogel.py index b5ab00a5..3b97addf 100644 --- a/samples/build_hydrogel.py +++ b/samples/build_hydrogel.py @@ -15,94 +15,241 @@ # # You should have received a copy of the GNU General Public License # along with this program. If not, see . +# +# Demonstrates: +# 1. Original hydrogel build (mpc=40, two bead types, no angles) with a 3D plot. +# 2. Angle-potential hydrogel (mpc=4, three bead types). +# Run with --include to also add crosslinker-adjacent angle templates. -import pyMBE +import argparse +import numpy as np import espressomd import matplotlib.pyplot as plt + +import pyMBE from pyMBE.lib.lattice import DiamondLattice -pmb = pyMBE.pymbe_library(seed=42) -reduced_unit_set = pmb.get_reduced_units() -# Monomers per chain -mpc = 40 -# Define node particle -NodeType = "node_type" -pmb.define_particle(name=NodeType, - sigma=0.355*pmb.units.nm, - epsilon=1*pmb.units('reduced_energy')) -# define monomers -BeadType1 = "C" -pmb.define_particle(name=BeadType1, - sigma=0.355*pmb.units.nm, - epsilon=1*pmb.units('reduced_energy')) -BeadType2 = "M" -pmb.define_particle(name=BeadType2, - sigma=0.355*pmb.units.nm, - epsilon=1*pmb.units('reduced_energy')) - -Res1 = "res_1" -pmb.define_residue(name=Res1, # Name of the residue - central_bead=BeadType1, # Define the central bead name - side_chains=[]) # Assuming no side chains for the monomer - - -Res2 = "res_2" -pmb.define_residue(name=Res2, # Name of the residue - central_bead=BeadType2, # Define the central bead name - side_chains=[]) # Assuming no side chains for the monomer - - -residue_list = [Res1]*(mpc//2) + [Res2]*(mpc//2) -pmb.define_molecule(name="hydrogel_chain", - residue_list=residue_list) - - -# Defining bonds in the hydrogel for all different pairs -generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2') -generic_bond_l = 0.355*pmb.units.nm -HARMONIC_parameters = {'r_0' : generic_bond_l, - 'k' : generic_harmonic_constant} -pmb.define_bond(bond_type = 'harmonic', - bond_parameters = HARMONIC_parameters, particle_pairs = [[BeadType1, BeadType1], - [BeadType1, BeadType2], - [BeadType2, BeadType2], - [NodeType, BeadType1], - [NodeType, BeadType2]]) -# Provide mpc and bond_l to Diamond Lattice -diamond_lattice = DiamondLattice(mpc, generic_bond_l) -espresso_system = espressomd.System(box_l = [diamond_lattice.box_l]*3) - -lattice_builder = pmb.initialize_lattice_builder(diamond_lattice) - -# Setting up node topology -indices = diamond_lattice.indices -node_topology = [] - -for index in range(len(indices)): - node_topology.append({"particle_name": NodeType, - "lattice_index": indices[index]}) - -# Setting up chain topology -connectivity = diamond_lattice.connectivity -node_labels = lattice_builder.node_labels -reverse_node_labels = {v: k for k, v in node_labels.items()} -connectivity_with_labels = {(reverse_node_labels[i], reverse_node_labels[j]) for i, j in connectivity} -chain_topology = [] - -for node_s, node_e in connectivity_with_labels: - chain_topology.append({'node_start':node_s, - 'node_end': node_e, - 'molecule_name':"hydrogel_chain"}) - -lattice_builder.chains = chain_topology - -pmb.define_hydrogel("my_hydrogel",node_topology, chain_topology) -hydrogel_info = pmb.create_hydrogel("my_hydrogel", espresso_system) - -fig = plt.figure() -ax = fig.add_subplot(111,projection="3d") -lattice_builder.draw_lattice(ax=ax, - pmb=pmb) -lattice_builder.draw_simulation_box(ax) -plt.legend(fontsize=12) -plt.show() + +def run_original_hydrogel(espresso_system): + """Build and visualize the original hydrogel (mpc=40, no angles).""" + pmb = pyMBE.pymbe_library(seed=42) + mpc = 40 + NodeType = "node_type" + pmb.define_particle(name=NodeType, + sigma=0.355*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + BeadType1 = "C" + pmb.define_particle(name=BeadType1, + sigma=0.355*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + BeadType2 = "M" + pmb.define_particle(name=BeadType2, + sigma=0.355*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + + Res1 = "res_1" + pmb.define_residue(name=Res1, central_bead=BeadType1, side_chains=[]) + Res2 = "res_2" + pmb.define_residue(name=Res2, central_bead=BeadType2, side_chains=[]) + + residue_list = [Res1]*(mpc//2) + [Res2]*(mpc//2) + pmb.define_molecule(name="hydrogel_chain", residue_list=residue_list) + + generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2') + generic_bond_l = 0.355*pmb.units.nm + HARMONIC_parameters = {'r_0': generic_bond_l, 'k': generic_harmonic_constant} + pmb.define_bond(bond_type='harmonic', + bond_parameters=HARMONIC_parameters, + particle_pairs=[[BeadType1, BeadType1], + [BeadType1, BeadType2], + [BeadType2, BeadType2], + [NodeType, BeadType1], + [NodeType, BeadType2]]) + + diamond_lattice = DiamondLattice(mpc, generic_bond_l) + espresso_system.box_l = [diamond_lattice.box_l] * 3 + lattice_builder = pmb.initialize_lattice_builder(diamond_lattice) + + indices = diamond_lattice.indices + node_topology = [] + for index in range(len(indices)): + node_topology.append({"particle_name": NodeType, + "lattice_index": indices[index]}) + + connectivity = diamond_lattice.connectivity + node_labels = lattice_builder.node_labels + reverse_node_labels = {v: k for k, v in node_labels.items()} + connectivity_with_labels = {(reverse_node_labels[i], reverse_node_labels[j]) + for i, j in connectivity} + chain_topology = [] + for node_s, node_e in connectivity_with_labels: + chain_topology.append({'node_start': node_s, + 'node_end': node_e, + 'molecule_name': "hydrogel_chain"}) + + lattice_builder.chains = chain_topology + pmb.define_hydrogel("my_hydrogel", node_topology, chain_topology) + pmb.create_hydrogel("my_hydrogel", espresso_system) + + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + lattice_builder.draw_lattice(ax=ax, pmb=pmb) + lattice_builder.draw_simulation_box(ax) + plt.legend(fontsize=12) + plt.show() + + +def define_angle_hydrogel_system(espresso_system, include_crosslinker_angles): + """Build the angle-potential hydrogel (mpc=4).""" + pmb = pyMBE.pymbe_library(seed=42) + + node_type = "D" + internal_bead_type = "A" + chain_end_bead_type = "C" + terminal_residue_name = "TerminalRes" + internal_residue_name = "InternalRes" + molecule_name = "hydrogel_chain" + + bond_length = 0.355 * pmb.units.nm + mpc = 4 + + pmb.define_particle(name=node_type, + sigma=bond_length, + epsilon=1*pmb.units("reduced_energy")) + pmb.define_particle(name=internal_bead_type, + sigma=bond_length, + epsilon=1*pmb.units("reduced_energy")) + pmb.define_particle(name=chain_end_bead_type, + sigma=bond_length, + epsilon=1*pmb.units("reduced_energy")) + + pmb.define_residue(name=terminal_residue_name, + central_bead=chain_end_bead_type, + side_chains=[]) + pmb.define_residue(name=internal_residue_name, + central_bead=internal_bead_type, + side_chains=[]) + pmb.define_molecule(name=molecule_name, + residue_list=[ + terminal_residue_name, + internal_residue_name, + internal_residue_name, + terminal_residue_name, + ]) + + harmonic_bond_parameters = { + "r_0": bond_length, + "k": 400 * pmb.units("reduced_energy / reduced_length**2"), + } + pmb.define_bond(bond_type="harmonic", + bond_parameters=harmonic_bond_parameters, + particle_pairs=[ + [internal_bead_type, internal_bead_type], + [internal_bead_type, chain_end_bead_type], + [node_type, chain_end_bead_type], + ]) + + angle_parameters = { + "k": 25 * pmb.units("reduced_energy"), + "phi_0": np.pi * pmb.units(""), + } + pmb.define_angle(angle_type="harmonic", + angle_parameters=angle_parameters, + particle_triplets=[(chain_end_bead_type, internal_bead_type, internal_bead_type)]) + + if include_crosslinker_angles: + pmb.define_angle(angle_type="harmonic", + angle_parameters=angle_parameters, + particle_triplets=[ + (node_type, chain_end_bead_type, internal_bead_type), # canonical: A-C-D + (chain_end_bead_type, node_type, chain_end_bead_type), # canonical: C-D-C + ]) + + diamond_lattice = DiamondLattice(mpc, bond_length) + espresso_system.box_l = [diamond_lattice.box_l] * 3 + pmb.initialize_lattice_builder(diamond_lattice) + + node_topology = [ + {"particle_name": node_type, "lattice_index": index} + for index in diamond_lattice.indices + ] + chain_topology = [] + for node_connectivity in diamond_lattice.connectivity: + node_start = str(diamond_lattice.indices[node_connectivity[0]]) + node_end = str(diamond_lattice.indices[node_connectivity[1]]) + chain_topology.append({ + "node_start": node_start, + "node_end": node_end, + "molecule_name": molecule_name, + }) + + pmb.define_hydrogel(name="simple_hydrogel", + node_map=node_topology, + chain_map=chain_topology) + hydrogel_id = pmb.create_hydrogel(name="simple_hydrogel", + espresso_system=espresso_system, + gen_angle=True) + return pmb, hydrogel_id + + +def summarize_angles(pmb): + """Return angle templates, instances, and counts by canonical angle name.""" + angle_templates_df = pmb.get_templates_df(pmb_type="angle") + angle_instances_df = pmb.get_instances_df(pmb_type="angle") + angle_counts = angle_instances_df["name"].value_counts().sort_index() + return angle_templates_df, angle_instances_df, angle_counts + + +def run_angle_case(label, include_crosslinker_angles, expected_angle_names, espresso_system): + """Clear the system, build one angle-hydrogel case, and assert expected angle names.""" + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + espresso_system.non_bonded_inter.reset() + + print(f"\n############ {label} ############") + pmb, hydrogel_id = define_angle_hydrogel_system(espresso_system, include_crosslinker_angles) + angle_templates_df, angle_instances_df, angle_counts = summarize_angles(pmb) + print(f"Hydrogel assembly id: {hydrogel_id}") + print("\nAngle templates:") + print(angle_templates_df) + print("\nAngle instance dataframe:") + print(angle_instances_df) + print("Angle instance counts:") + print(angle_counts) + + angle_names = set(angle_counts.index) + if angle_names != expected_angle_names: + raise AssertionError( + f"{label}: expected angle names {sorted(expected_angle_names)}, " + f"got {sorted(angle_names)}" + ) + print(f"{label}: OK") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Build hydrogels: original (mpc=40) then angle-potential (mpc=4)." + ) + parser.add_argument( + "--include", + action="store_true", + help="Include crosslinker-adjacent angles (D-C-A and C-D-C) in the angle demo.", + ) + args = parser.parse_args() + + espresso_system = espressomd.System(box_l=[1.0] * 3) + + # Part 1: original hydrogel (mpc=40, no angles) + run_original_hydrogel(espresso_system) + + # Part 2: angle-potential hydrogel (mpc=4); clearing is done inside run_angle_case + if args.include: + run_angle_case(label="Chain + crosslinker angles", + include_crosslinker_angles=True, + expected_angle_names={"A-A-C", "A-C-D", "C-D-C"}, + espresso_system=espresso_system) + else: + run_angle_case(label="Chain angles only", + include_crosslinker_angles=False, + expected_angle_names={"A-A-C"}, + espresso_system=espresso_system) diff --git a/samples/build_hydrogel_angle.py b/samples/build_hydrogel_angle.py deleted file mode 100644 index bae855ec..00000000 --- a/samples/build_hydrogel_angle.py +++ /dev/null @@ -1,161 +0,0 @@ -# -# Example: demonstrating hydrogel angle generation in pyMBE. -# -# This script builds a simple hydrogel of: -# - N: crosslinker/node particles -# - C: residue/chain particles -# -# By default, it defines only C-C-C and generates only chain-internal angles. -# If run with `--include`, it also defines the crosslinker-adjacent -# angle templates N-C-C and C-N-C. -# -# Note: -# - The canonical angle name for N-C-C is C-C-N, because pyMBE sorts the -# two side particles alphabetically when building angle template names. -# - -import argparse -import numpy as np -import espressomd - -import pyMBE -from pyMBE.lib.lattice import DiamondLattice - - -def define_simple_hydrogel_system(include_crosslinker_angles): - """ - Build a simple hydrogel with C residues and N crosslinkers. - """ - pmb = pyMBE.pymbe_library(seed=42) - - node_type = "N" - bead_type = "C" - residue_name = "Res" - molecule_name = "hydrogel_chain" - - bond_length = 0.355 * pmb.units.nm - mpc = 4 - - pmb.define_particle(name=node_type, - sigma=bond_length, - epsilon=1 * pmb.units("reduced_energy")) - pmb.define_particle(name=bead_type, - sigma=bond_length, - epsilon=1 * pmb.units("reduced_energy")) - - pmb.define_residue(name=residue_name, - central_bead=bead_type, - side_chains=[]) - pmb.define_molecule(name=molecule_name, - residue_list=[residue_name] * mpc) - - harmonic_bond_parameters = { - "r_0": bond_length, - "k": 400 * pmb.units("reduced_energy / reduced_length**2"), - } - pmb.define_bond(bond_type="harmonic", - bond_parameters=harmonic_bond_parameters, - particle_pairs=[ - [bead_type, bead_type], - [node_type, bead_type], - ]) - - angle_parameters = { - "k": 25 * pmb.units("reduced_energy"), - "phi_0": np.pi * pmb.units(""), - } - pmb.define_angle(angle_type="harmonic", - angle_parameters=angle_parameters, - particle_triplets=[(bead_type, bead_type, bead_type)]) - - if include_crosslinker_angles: - pmb.define_angle(angle_type="harmonic", - angle_parameters=angle_parameters, - particle_triplets=[ - (node_type, bead_type, bead_type), # canonical: C-C-N - (bead_type, node_type, bead_type), # canonical: C-N-C - ]) - - diamond_lattice = DiamondLattice(mpc, bond_length) - espresso_system = espressomd.System(box_l=[diamond_lattice.box_l] * 3) - pmb.initialize_lattice_builder(diamond_lattice) - - node_topology = [ - {"particle_name": node_type, "lattice_index": index} - for index in diamond_lattice.indices - ] - chain_topology = [] - for node_connectivity in diamond_lattice.connectivity: - node_start = str(diamond_lattice.indices[node_connectivity[0]]) - node_end = str(diamond_lattice.indices[node_connectivity[1]]) - chain_topology.append({ - "node_start": node_start, - "node_end": node_end, - "molecule_name": molecule_name, - }) - - pmb.define_hydrogel(name="simple_hydrogel", - node_map=node_topology, - chain_map=chain_topology) - hydrogel_id = pmb.create_hydrogel(name="simple_hydrogel", - espresso_system=espresso_system, - gen_angle=True) - return pmb, espresso_system, hydrogel_id - - -def summarize_angles(pmb): - """ - Return angle templates, angle instances, and counts by canonical angle name. - """ - angle_templates_df = pmb.get_templates_df(pmb_type="angle") - angle_instances_df = pmb.get_instances_df(pmb_type="angle") - angle_counts = angle_instances_df["name"].value_counts().sort_index() - return angle_templates_df, angle_instances_df, angle_counts - - -def run_case(label, include_crosslinker_angles, expected_angle_names): - """ - Build one hydrogel case, print angle counts, and assert expectations. - """ - print(f"\n############ {label} ############") - pmb, _, hydrogel_id = define_simple_hydrogel_system( - include_crosslinker_angles=include_crosslinker_angles - ) - angle_templates_df, angle_instances_df, angle_counts = summarize_angles(pmb) - print(f"Hydrogel assembly id: {hydrogel_id}") - print("\nAngle templates:") - print(angle_templates_df) - print("\nAngle instance dataframe:") - print(angle_instances_df) - print("Angle instance counts:") - print(angle_counts) - - angle_names = set(angle_counts.index) - if angle_names != expected_angle_names: - raise AssertionError( - f"{label}: expected angle names {sorted(expected_angle_names)}, " - f"got {sorted(angle_names)}" - ) - - print(f"{label}: OK") - - -if __name__ == "__main__": - parser = argparse.ArgumentParser( - description="Build a simple hydrogel and optionally include crosslinker-adjacent angles." - ) - parser.add_argument( - "--include", - action="store_true", - help="Include crosslinker-adjacent angles (N-C-C and C-N-C).", - ) - args = parser.parse_args() - - if args.include: - run_case(label="Chain + crosslinker angles", - include_crosslinker_angles=True, - expected_angle_names={"C-C-C", "C-C-N", "C-N-C"}) - else: - run_case(label="Chain angles only", - include_crosslinker_angles=False, - expected_angle_names={"C-C-C"}) diff --git a/testsuite/angle_tests.py b/testsuite/angle_tests.py index 68a6afa8..072dd9be 100644 --- a/testsuite/angle_tests.py +++ b/testsuite/angle_tests.py @@ -1,5 +1,5 @@ # -# Copyright (C) 2024-2026 pyMBE-dev team +# Copyright (C) 2026 pyMBE-dev team # # This file is part of pyMBE. # diff --git a/testsuite/hydrogel_builder_with_angles.py b/testsuite/hydrogel_builder_with_angles.py index 658c89d5..4f91c62e 100644 --- a/testsuite/hydrogel_builder_with_angles.py +++ b/testsuite/hydrogel_builder_with_angles.py @@ -1,5 +1,5 @@ # -# Copyright (C) 2024-2026 pyMBE-dev team +# Copyright (C) 2026 pyMBE-dev team # # This file is part of pyMBE. # From 059a8c3233a7f13e8f17958b7ed40631795535cc Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sun, 3 May 2026 18:56:11 +0200 Subject: [PATCH 09/13] doc strings improved --- pyMBE/pyMBE.py | 10 +- samples/build_hydrogel.py | 90 ++++++++-- samples/setup_angular_potential.py | 200 ++++++++++++++++++++++ testsuite/angle_tests.py | 58 +++---- testsuite/hydrogel_builder_with_angles.py | 51 +++++- testsuite/test_io_database.py | 6 +- 6 files changed, 359 insertions(+), 56 deletions(-) create mode 100644 samples/setup_angular_potential.py diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index 99a7098e..22746bc0 100644 --- a/pyMBE/pyMBE.py +++ b/pyMBE/pyMBE.py @@ -372,7 +372,7 @@ def _generate_hydrogel_crosslinker_angles(self, espresso_system, central_particl ) for side_particle_id1, central_particle_id, side_particle_id3, _ in triplets: - self.create_angle(particle_id1=side_particle_id1, + self.create_angular_potential(particle_id1=side_particle_id1, particle_id2=central_particle_id, particle_id3=side_particle_id3, espresso_system=espresso_system, @@ -1575,7 +1575,7 @@ def define_default_bond(self, bond_type, bond_parameters): tpl.name = "default" self.db._register_template(tpl) - def define_angle(self, angle_type, angle_parameters, particle_triplets): + def define_angular_potential(self, angle_type, angle_parameters, particle_triplets): """ Defines angle potential templates for each particle triplet in `particle_triplets`. @@ -1622,7 +1622,7 @@ def define_angle(self, angle_type, angle_parameters, particle_triplets): angle_names.append(tpl.name) self.db._register_template(tpl) - def define_default_angle(self, angle_type, angle_parameters): + def define_default_angular_potential(self, angle_type, angle_parameters): """ Defines an angle template as a "default" template in the pyMBE database. @@ -1655,7 +1655,7 @@ def define_default_angle(self, angle_type, angle_parameters): tpl.name = "default" self.db._register_template(tpl) - def create_angle(self, particle_id1, particle_id2, particle_id3, espresso_system, use_default_angle=False): + def create_angular_potential(self, particle_id1, particle_id2, particle_id3, espresso_system, use_default_angle=False): """ Creates an angle between three particle instances in an ESPResSo system and registers it in the pyMBE database. @@ -1805,7 +1805,7 @@ def _generate_angles_for_entity(self, espresso_system, entity_id, entity_id_col) i = nbs[idx_i] k = nbs[idx_k] try: - self.create_angle(particle_id1=i, + self.create_angular_potential(particle_id1=i, particle_id2=j, particle_id3=k, espresso_system=espresso_system, diff --git a/samples/build_hydrogel.py b/samples/build_hydrogel.py index 3b97addf..56f430b5 100644 --- a/samples/build_hydrogel.py +++ b/samples/build_hydrogel.py @@ -30,8 +30,17 @@ from pyMBE.lib.lattice import DiamondLattice -def run_original_hydrogel(espresso_system): - """Build and visualize the original hydrogel (mpc=40, no angles).""" +def run_simple_hydrogel(espresso_system): + """ + Build and visualize a simple hydrogel with no angular potential (mpc=40). + + Defines two bead types (C, M) and one node type, builds the diamond-lattice + hydrogel, and displays a 3D plot of the lattice. + + Args: + espresso_system (espressomd.system.System): ESPResSo system object in + which the hydrogel particles and bonds will be created. + """ pmb = pyMBE.pymbe_library(seed=42) mpc = 40 NodeType = "node_type" @@ -99,8 +108,29 @@ def run_original_hydrogel(espresso_system): plt.show() -def define_angle_hydrogel_system(espresso_system, include_crosslinker_angles): - """Build the angle-potential hydrogel (mpc=4).""" +def define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_angles): + """ + Build a hydrogel with angular potentials on the diamond lattice (mpc=4). + + Defines three bead types: crosslinker nodes (D), chain-end beads (C), and + internal chain beads (A). Always defines the chain-internal angle A-A-C. + When include_crosslinker_angles is True, also defines D-C-A and C-D-C angles. + + Args: + espresso_system (espressomd.system.System): ESPResSo system object in + which the hydrogel particles, bonds, and angles will be created. + Its box_l is updated to match the diamond lattice dimensions. + include_crosslinker_angles (bool): If True, also define angle templates + for the crosslinker-adjacent triplets D-C-A (canonical: A-C-D) and + C-D-C in addition to the chain-internal A-A-C angle. + + Returns: + tuple: + - pmb (pyMBE.pymbe_library): pyMBE instance holding all templates + and instances for this hydrogel. + - hydrogel_id (int): Assembly ID of the created hydrogel, as + returned by pmb.create_hydrogel. + """ pmb = pyMBE.pymbe_library(seed=42) node_type = "D" @@ -153,12 +183,12 @@ def define_angle_hydrogel_system(espresso_system, include_crosslinker_angles): "k": 25 * pmb.units("reduced_energy"), "phi_0": np.pi * pmb.units(""), } - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=angle_parameters, particle_triplets=[(chain_end_bead_type, internal_bead_type, internal_bead_type)]) if include_crosslinker_angles: - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=angle_parameters, particle_triplets=[ (node_type, chain_end_bead_type, internal_bead_type), # canonical: A-C-D @@ -193,21 +223,51 @@ def define_angle_hydrogel_system(espresso_system, include_crosslinker_angles): def summarize_angles(pmb): - """Return angle templates, instances, and counts by canonical angle name.""" + """ + Return angle templates, instances, and per-name counts from a pyMBE instance. + + Args: + pmb (pyMBE.pymbe_library): pyMBE instance whose database will be queried. + + Returns: + tuple: + - angle_templates_df (pandas.DataFrame): DataFrame of all defined + angle templates. + - angle_instances_df (pandas.DataFrame): DataFrame of all angle + instances currently in the system. + - angle_counts (pandas.Series): Count of angle instances grouped + by canonical angle name, sorted alphabetically by name. + """ angle_templates_df = pmb.get_templates_df(pmb_type="angle") angle_instances_df = pmb.get_instances_df(pmb_type="angle") angle_counts = angle_instances_df["name"].value_counts().sort_index() return angle_templates_df, angle_instances_df, angle_counts -def run_angle_case(label, include_crosslinker_angles, expected_angle_names, espresso_system): - """Clear the system, build one angle-hydrogel case, and assert expected angle names.""" +def run_hydrogel_with_angular_potential(label, include_crosslinker_angles, expected_angle_names, espresso_system): + """ + Clear the ESPResSo system, build one angle-hydrogel case, and assert the + expected set of canonical angle names is generated. + + Args: + label (str): Human-readable name for this case, printed as a section + header and in assertion error messages. + include_crosslinker_angles (bool): Passed directly to + define_hydrogel_with_angular_potential; controls whether + crosslinker-adjacent angle templates are included. + expected_angle_names (set of str): Set of canonical angle names + (e.g. ``{"A-A-C", "A-C-D"}``}) that must be present in the angle + instance DataFrame. An AssertionError is raised if the actual set + differs. + espresso_system (espressomd.system.System): ESPResSo system object to + clear and reuse for this case. + """ espresso_system.part.clear() espresso_system.bonded_inter.clear() espresso_system.non_bonded_inter.reset() print(f"\n############ {label} ############") - pmb, hydrogel_id = define_angle_hydrogel_system(espresso_system, include_crosslinker_angles) + pmb, hydrogel_id = define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_angles) angle_templates_df, angle_instances_df, angle_counts = summarize_angles(pmb) print(f"Hydrogel assembly id: {hydrogel_id}") print("\nAngle templates:") @@ -239,17 +299,17 @@ def run_angle_case(label, include_crosslinker_angles, expected_angle_names, espr espresso_system = espressomd.System(box_l=[1.0] * 3) - # Part 1: original hydrogel (mpc=40, no angles) - run_original_hydrogel(espresso_system) + # Part 1: simple hydrogel (mpc=40, no angles) + run_simple_hydrogel(espresso_system) - # Part 2: angle-potential hydrogel (mpc=4); clearing is done inside run_angle_case + # Part 2: angle-potential hydrogel (mpc=4); clearing is done inside run_hydrogel_with_angular_potential if args.include: - run_angle_case(label="Chain + crosslinker angles", + run_hydrogel_with_angular_potential(label="Chain + crosslinker angles", include_crosslinker_angles=True, expected_angle_names={"A-A-C", "A-C-D", "C-D-C"}, espresso_system=espresso_system) else: - run_angle_case(label="Chain angles only", + run_hydrogel_with_angular_potential(label="Chain angles only", include_crosslinker_angles=False, expected_angle_names={"A-A-C"}, espresso_system=espresso_system) diff --git a/samples/setup_angular_potential.py b/samples/setup_angular_potential.py new file mode 100644 index 00000000..80183fab --- /dev/null +++ b/samples/setup_angular_potential.py @@ -0,0 +1,200 @@ +# +# 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 . +# +# Example: demonstrating angle potential support in pyMBE. +# +# This script shows three usages of `gen_angle=True`: +# 1. create_residue with gen_angle=True +# - builds a single residue with auto-generated angles +# 2. create_molecule with gen_angle=True +# - builds a chain of residues; angles are auto-generated across the +# full molecule, including angles that span residue boundaries +# 3. create_residue with a nested residue as side chain +# - builds a residue whose side chain is itself a residue; verifies +# that all bonded triplets (intra- and cross-level) are generated + +import pyMBE +import numpy as np +import espressomd + +# ---------------------------------------------------------------------- +# Set up the ESPResSo system and pyMBE +# ---------------------------------------------------------------------- +espresso_system = espressomd.System(box_l=[20] * 3) +pmb = pyMBE.pymbe_library(seed=42) + +# ---------------------------------------------------------------------- +# Particle definitions +# ---------------------------------------------------------------------- +pmb.define_particle(name='A', + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + +pmb.define_particle(name='B', + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + +pmb.define_particle(name='C', + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + +pmb.define_particle(name='D', + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + +# ---------------------------------------------------------------------- +# Default bond used everywhere +# ---------------------------------------------------------------------- +HARMONIC_bond_parameters = { + 'r_0': 0.4 * pmb.units.nm, + 'k': 40 * pmb.units('reduced_energy / reduced_length**2'), +} +pmb.define_default_bond(bond_type='harmonic', + bond_parameters=HARMONIC_bond_parameters) + +# ---------------------------------------------------------------------- +# Default angle used by both demos +# ---------------------------------------------------------------------- +default_angle_params = { + 'k': 50 * pmb.units('reduced_energy'), + 'phi_0': np.pi * pmb.units(''), +} +pmb.define_default_angular_potential(angle_type="harmonic", + angle_parameters=default_angle_params) + + +def show_database(label): + print(f"\n############ {label} ############") + for pmb_type in ["particle", "residue", "molecule", "bond", "angle"]: + print(f"\n=== {pmb_type} templates ===") + print(pmb.get_templates_df(pmb_type=pmb_type)) + print(f"\n=== {pmb_type} instances ===") + print(pmb.get_instances_df(pmb_type=pmb_type)) + + +# ====================================================================== +# Demo 1: create_residue with gen_angle=True +# ====================================================================== +# Topology: central A bonded to side particles A, B, C. +# Auto-generated angles (central A has 3 neighbors): A-A-B, A-A-C, B-A-C +pmb.define_residue(name="Res_1", + central_bead="A", + side_chains=["A", "B", "C"]) + +residue_id = pmb.create_residue(name="Res_1", + espresso_system=espresso_system, + use_default_bond=True, + gen_angle=True) + +show_database("Demo 1: single residue with gen_angle=True") + +# ---------------------------------------------------------------------- +# Clean up: remove the residue (and its particles, bonds, angles) from +# both ESPResSo and the pyMBE database before running Demo 2. +# ---------------------------------------------------------------------- +pmb.delete_instances_in_system(instance_id=residue_id, + pmb_type="residue", + espresso_system=espresso_system) + + +# ====================================================================== +# Demo 2: create_residue with a nested residue as side chain +# ====================================================================== +# Topology: D₁ where SubRes_3 = B bonded to C and D₂ +# \ +# A +# | +# B -- C +# | +# D₂ +# +# SubRes_3 (central B, side chains C and D) is passed as a side chain of +# the outer residue NestedRes_3 (central A, direct side chain D). +# D₁ and D₂ are distinct instances of particle type D. +# +# All bonded triplets and their canonical angle names: +# centered at A: D₁ -- A -- B → "B-A-D" (outer {B,D} sorted) +# centered at B: A -- B -- C → "A-B-C" +# centered at B: A -- B -- D₂ → "A-B-D" +# centered at B: C -- B -- D₂ → "C-B-D" +pmb.define_residue(name="SubRes_3", + central_bead="B", + side_chains=["C", "D"]) + +pmb.define_residue(name="NestedRes_3", + central_bead="A", + side_chains=["D", "SubRes_3"]) + +nested_residue_id = pmb.create_residue(name="NestedRes_3", + espresso_system=espresso_system, + use_default_bond=True, + gen_angle=True) + +show_database("Demo 2: nested residue with gen_angle=True") + +angle_instances_df = pmb.get_instances_df(pmb_type="angle") +particle_instances_df = pmb.get_instances_df(pmb_type="particle") +print(f"\nAngle instance count: {len(angle_instances_df)}") + +id_to_type = dict(zip(particle_instances_df["particle_id"], + particle_instances_df["name"])) +actual_triplets = set() +for _, row in angle_instances_df.iterrows(): + p1 = id_to_type[row["particle_id1"]] + center = id_to_type[row["particle_id2"]] + p3 = id_to_type[row["particle_id3"]] + outer = sorted([p1, p3]) + actual_triplets.add(f"{outer[0]}-{center}-{outer[1]}") + +print("Resolved angle triplets:", sorted(actual_triplets)) +expected_triplets = {"A-B-C", "A-B-D", "B-A-D", "C-B-D"} +assert actual_triplets == expected_triplets, ( + f"Demo 2: expected {sorted(expected_triplets)}, got {sorted(actual_triplets)}" +) +print("Demo 2: all expected angle triplets generated OK") + +# ====================================================================== +# Demo 3: create_molecule with gen_angle=True +# ====================================================================== +# Topology: a chain of 4 Res_2 residues, each with central A and one +# side chain B. After auto-generation, angles include both intra-residue +# (A-A-B) and inter-residue (A-A-A) triplets. +pmb.define_residue(name="Res_2", + central_bead="A", + side_chains=["B"]) + +pmb.define_molecule(name="Mol_1", + residue_list=["Res_2"] * 4) + +molecule_ids = pmb.create_molecule(name="Mol_1", + number_of_molecules=1, + espresso_system=espresso_system, + backbone_vector=[1.0, 0.0, 0.0], + use_default_bond=True, + gen_angle=True) + +show_database("Demo 3: linear molecule with gen_angle=True") + +pmb.delete_instances_in_system(instance_id=molecule_ids[0], + pmb_type="molecule", + espresso_system=espresso_system) + diff --git a/testsuite/angle_tests.py b/testsuite/angle_tests.py index 072dd9be..14f97f8f 100644 --- a/testsuite/angle_tests.py +++ b/testsuite/angle_tests.py @@ -119,7 +119,7 @@ def test_angle_setup_harmonic(self): particle_pairs=[['A', 'B'], ['B', 'C']]) # Define harmonic angle for triplet A-B-C - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('A', 'B', 'C')]) @@ -143,7 +143,7 @@ def test_angle_setup_harmonic(self): espresso_system=espresso_system) # Create the angle: A-B-C (B is central) - pmb.create_angle(particle_id1=pid_A[0], + pmb.create_angular_potential(particle_id1=pid_A[0], particle_id2=pid_B[0], particle_id3=pid_C[0], espresso_system=espresso_system) @@ -176,7 +176,7 @@ def test_angle_setup_cosine(self): particle_pairs=[['A', 'B']]) # Define cosine angle for triplet A-B-A - pmb.define_angle(angle_type="cosine", + pmb.define_angular_potential(angle_type="cosine", angle_parameters=self.cosine_angle_params, particle_triplets=[('A', 'B', 'A')]) @@ -199,7 +199,7 @@ def test_angle_setup_cosine(self): particle_id2=pid_A2[0], espresso_system=espresso_system) - pmb.create_angle(particle_id1=pid_A1[0], + pmb.create_angular_potential(particle_id1=pid_A1[0], particle_id2=pid_B[0], particle_id3=pid_A2[0], espresso_system=espresso_system) @@ -230,7 +230,7 @@ def test_angle_setup_harmonic_cosine(self): bond_parameters=self.harmonic_bond_params, particle_pairs=[['A', 'B'], ['B', 'C']]) - pmb.define_angle(angle_type="harmonic_cosine", + pmb.define_angular_potential(angle_type="harmonic_cosine", angle_parameters=self.harmonic_cosine_angle_params, particle_triplets=[('A', 'B', 'C')]) @@ -251,7 +251,7 @@ def test_angle_setup_harmonic_cosine(self): particle_id2=pid_C[0], espresso_system=espresso_system) - pmb.create_angle(particle_id1=pid_A[0], + pmb.create_angular_potential(particle_id1=pid_A[0], particle_id2=pid_B[0], particle_id3=pid_C[0], espresso_system=espresso_system) @@ -277,7 +277,7 @@ def test_get_espresso_angle_instance_reuses_cached_object(self): pmb = pyMBE.pymbe_library(seed=42) self.define_templates(pmb) - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('A', 'B', 'C')]) @@ -308,7 +308,7 @@ def test_default_angle(self): particle_pairs=[['A', 'B'], ['B', 'C']]) # Define only a default angle (no specific triplet) - pmb.define_default_angle(angle_type="harmonic", + pmb.define_default_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params) # Create three particles @@ -331,7 +331,7 @@ def test_default_angle(self): espresso_system=espresso_system) # Create angle using default (no specific A-B-C template exists) - pmb.create_angle(particle_id1=pid_A[0], + pmb.create_angular_potential(particle_id1=pid_A[0], particle_id2=pid_B[0], particle_id3=pid_C[0], espresso_system=espresso_system, @@ -365,11 +365,11 @@ def test_specific_angle_over_default(self): particle_pairs=[['A', 'B'], ['B', 'C']]) # Define default angle with cosine parameters - pmb.define_default_angle(angle_type="cosine", + pmb.define_default_angular_potential(angle_type="cosine", angle_parameters=self.cosine_angle_params) # Define specific angle for A-B-C with harmonic parameters - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('A', 'B', 'C')]) @@ -392,7 +392,7 @@ def test_specific_angle_over_default(self): espresso_system=espresso_system) # Should use the specific A-B-C template, not the default - pmb.create_angle(particle_id1=pid_A[0], + pmb.create_angular_potential(particle_id1=pid_A[0], particle_id2=pid_B[0], particle_id3=pid_C[0], espresso_system=espresso_system, @@ -421,7 +421,7 @@ def test_angle_raised_exceptions(self): pmb = pyMBE.pymbe_library(seed=42) self.define_templates(pmb) - for callback in [pmb.define_angle, pmb.define_default_angle]: + for callback in [pmb.define_angular_potential, pmb.define_default_angular_potential]: with self.subTest(msg=f'using method {callback.__qualname__}()'): self.check_angle_exceptions(callback, pmb) @@ -432,7 +432,7 @@ def check_angle_exceptions(self, callback, pmb): 'phi_0': np.pi * pmb.units('')} input_parameters = {"angle_type": angle_type, "angle_parameters": angle_params} - if callback == pmb.define_angle: + if callback == pmb.define_angular_potential: input_parameters["particle_triplets"] = [('A', 'B', 'C')] np.testing.assert_raises(NotImplementedError, callback, **input_parameters) @@ -442,7 +442,7 @@ def check_angle_exceptions(self, callback, pmb): angle_params = {'phi_0': np.pi * pmb.units('')} input_parameters = {"angle_type": angle_type, "angle_parameters": angle_params} - if callback == pmb.define_angle: + if callback == pmb.define_angular_potential: input_parameters["particle_triplets"] = [('A', 'B', 'C')] np.testing.assert_raises(ValueError, callback, **input_parameters) @@ -452,19 +452,19 @@ def check_angle_exceptions(self, callback, pmb): angle_params = {'k': 50 * pmb.units('reduced_energy')} input_parameters = {"angle_type": angle_type, "angle_parameters": angle_params} - if callback == pmb.define_angle: + if callback == pmb.define_angular_potential: input_parameters["particle_triplets"] = [('A', 'B', 'C')] np.testing.assert_raises(ValueError, callback, **input_parameters) - # Check exception for duplicate triplet in define_angle - if callback == pmb.define_angle: + # Check exception for duplicate triplet in define_angular_potential + if callback == pmb.define_angular_potential: test = {"angle_type": "harmonic", "angle_parameters": self.harmonic_angle_params, "particle_triplets": [("A", "B", "A"), ("A", "B", "A")]} np.testing.assert_raises(RuntimeError, - pmb.define_angle, + pmb.define_angular_potential, **test) def test_sanity_get_angle_template(self): @@ -487,7 +487,7 @@ def test_canonical_angle_name(self): self.define_templates(pmb) # Define angle with C-B-A order — should produce same template as A-B-C - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('C', 'B', 'A')]) @@ -513,7 +513,7 @@ def test_angle_without_bonds_raises_error(self): self.define_templates(pmb) # Define angle template but no bond templates - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('A', 'B', 'C')]) @@ -529,8 +529,8 @@ def test_angle_without_bonds_raises_error(self): number_of_particles=1) # Attempting to create angle without bonds should raise ValueError - with self.assertRaises(ValueError, msg="create_angle should raise ValueError when no bonds exist"): - pmb.create_angle(particle_id1=pid_A[0], + with self.assertRaises(ValueError, msg="create_angular_potential should raise ValueError when no bonds exist"): + pmb.create_angular_potential(particle_id1=pid_A[0], particle_id2=pid_B[0], particle_id3=pid_C[0], espresso_system=espresso_system) @@ -554,7 +554,7 @@ def test_angle_with_partial_bonds_raises_error(self): bond_parameters=self.harmonic_bond_params, particle_pairs=[['A', 'B']]) - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('A', 'B', 'C')]) @@ -574,8 +574,8 @@ def test_angle_with_partial_bonds_raises_error(self): espresso_system=espresso_system) # Should raise ValueError because B-C bond is missing - with self.assertRaises(ValueError, msg="create_angle should raise ValueError when a bond is missing"): - pmb.create_angle(particle_id1=pid_A[0], + with self.assertRaises(ValueError, msg="create_angular_potential should raise ValueError when a bond is missing"): + pmb.create_angular_potential(particle_id1=pid_A[0], particle_id2=pid_B[0], particle_id3=pid_C[0], espresso_system=espresso_system) @@ -611,7 +611,7 @@ def test_generate_angles_for_entity_creates_angle_from_bonds(self): pmb.define_bond(bond_type="harmonic", bond_parameters=self.harmonic_bond_params, particle_pairs=[['A', 'B'], ['B', 'C']]) - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('A', 'B', 'C')]) @@ -710,7 +710,7 @@ def test_create_residue_with_gen_angle_generates_angles(self): pmb.define_bond(bond_type="harmonic", bond_parameters=self.harmonic_bond_params, particle_pairs=[['A', 'B'], ['B', 'C']]) - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('A', 'B', 'C')]) pmb.define_residue(name="R_angle", @@ -748,7 +748,7 @@ def test_create_molecule_with_gen_angle_generates_angles(self): pmb.define_bond(bond_type="harmonic", bond_parameters=self.harmonic_bond_params, particle_pairs=[['A', 'B'], ['B', 'C']]) - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=self.harmonic_angle_params, particle_triplets=[('A', 'B', 'C')]) pmb.define_residue(name="R_A", diff --git a/testsuite/hydrogel_builder_with_angles.py b/testsuite/hydrogel_builder_with_angles.py index 4f91c62e..632c1f74 100644 --- a/testsuite/hydrogel_builder_with_angles.py +++ b/testsuite/hydrogel_builder_with_angles.py @@ -29,6 +29,32 @@ def build_simple_hydrogel_with_optional_angles(junction_angle_mode="none", mpc_local=4): """ Build a simple hydrogel used to validate hydrogel-specific angle support. + + Defines two particle types: chain beads (C) and crosslinker nodes (N), + connected in a diamond lattice. Always defines the chain-internal C-C-C + angular potential template. Crosslinker-adjacent templates are added + depending on junction_angle_mode. + + Args: + junction_angle_mode (str): Controls which crosslinker angular potential + templates are defined. Accepted values: + - ``"none"``: no crosslinker angles defined (default). + - ``"partial"``: only C-N-C defined (missing N-C-C, triggers error + on hydrogel creation with gen_angle=True). + - ``"full"``: both C-N-C and C-C-N defined. + mpc_local (int): Number of monomers per chain in the diamond lattice. + Defaults to 4. + + Returns: + tuple: + - pmb_local (pyMBE.pymbe_library): pyMBE instance with all + templates and topology defined but hydrogel not yet created. + - espresso_system (espressomd.system.System): ESPResSo system + where the hydrogel will be created. + - hydrogel_name_local (str): Name under which the hydrogel is + registered in pmb_local (``"simple_hydrogel"``). + - diamond_lattice_local (DiamondLattice): The diamond lattice + instance used to build the topology. """ pmb_local = pyMBE.pymbe_library(seed=7) node_type = "N" @@ -61,7 +87,7 @@ def build_simple_hydrogel_with_optional_angles(junction_angle_mode="none", mpc_l "k": 5 * pmb_local.units("reduced_energy"), "phi_0": np.pi * pmb_local.units(""), } - pmb_local.define_angle(angle_type="harmonic", + pmb_local.define_angular_potential(angle_type="harmonic", angle_parameters=angle_parameters, particle_triplets=[(bead_type, bead_type, bead_type)]) @@ -69,7 +95,7 @@ def build_simple_hydrogel_with_optional_angles(junction_angle_mode="none", mpc_l particle_triplets = [(bead_type, node_type, bead_type)] if junction_angle_mode == "full": particle_triplets.append((node_type, bead_type, bead_type)) - pmb_local.define_angle(angle_type="harmonic", + pmb_local.define_angular_potential(angle_type="harmonic", angle_parameters=angle_parameters, particle_triplets=particle_triplets) @@ -100,7 +126,15 @@ def build_simple_hydrogel_with_optional_angles(junction_angle_mode="none", mpc_l def get_angle_counts(pmb_local): """ - Return counts of angle instances by canonical angle name. + Return counts of angular potential instances grouped by canonical angle name. + + Args: + pmb_local (pyMBE.pymbe_library): pyMBE instance whose database will + be queried for angular potential instances. + + Returns: + collections.Counter: Maps each canonical angle name (str) to the + number of instances of that angular potential in the system. """ return Counter( angle.name @@ -110,7 +144,16 @@ def get_angle_counts(pmb_local): def expected_crosslinker_angle_counts(diamond_lattice_local): """ - Return expected chain and crosslinker angle counts for the simple hydrogel. + Return the expected angular potential instance counts for the simple hydrogel + when all crosslinker templates are defined (junction_angle_mode="full"). + + Args: + diamond_lattice_local (DiamondLattice): The diamond lattice instance + used to build the hydrogel, providing mpc and connectivity info. + + Returns: + dict: Maps canonical angle name (str) to the expected instance count + (int). Keys are ``"C-C-C"``, ``"C-C-N"``, and ``"C-N-C"``. """ number_of_nodes = 8 number_of_chains = 16 diff --git a/testsuite/test_io_database.py b/testsuite/test_io_database.py index fcbd7b84..186eaa5c 100644 --- a/testsuite/test_io_database.py +++ b/testsuite/test_io_database.py @@ -303,11 +303,11 @@ def test_io_angle_templates(self): "phi_0": 1.0 * pmb.units("")} parameters2 = {"k": 30.0 * pmb.units.reduced_energy, "phi_0": 0.5 * pmb.units("")} - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=parameters1, particle_triplets=[("Z", "X", "Z"), ("X", "Z", "Y")]) - pmb.define_default_angle(angle_type="cosine", + pmb.define_default_angular_potential(angle_type="cosine", angle_parameters=parameters2) new_pmb = pyMBE.pymbe_library(23) with tempfile.TemporaryDirectory() as tmp_directory: @@ -489,7 +489,7 @@ def test_io_instances(self): ["X","X"]]) angle_parameters = {"k": 50.0 * pmb.units.reduced_energy, "phi_0": 1.0 * pmb.units("")} - pmb.define_angle(angle_type="harmonic", + pmb.define_angular_potential(angle_type="harmonic", angle_parameters=angle_parameters, particle_triplets=[("X", "Z", "Z")]) pmb.define_molecule(name="M1", From 2ca8b4454663cd8ebed04dabd45551ef361d9d02 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sun, 3 May 2026 19:09:21 +0200 Subject: [PATCH 10/13] rm build_angle.py --- samples/build_angle.py | 198 ----------------------------------------- 1 file changed, 198 deletions(-) delete mode 100644 samples/build_angle.py diff --git a/samples/build_angle.py b/samples/build_angle.py deleted file mode 100644 index c1465111..00000000 --- a/samples/build_angle.py +++ /dev/null @@ -1,198 +0,0 @@ -# -# 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 . -# -# Example: demonstrating angle potential support in pyMBE. -# -# This script shows three usages of `gen_angle=True`: -# 1. create_residue with gen_angle=True -# - builds a single residue with auto-generated angles -# 2. create_molecule with gen_angle=True -# - builds a chain of residues; angles are auto-generated across the -# full molecule, including angles that span residue boundaries -# 3. create_residue with a nested residue as side chain -# - builds a residue whose side chain is itself a residue; verifies -# that all bonded triplets (intra- and cross-level) are generated - -import pyMBE -import numpy as np -import espressomd - -# ---------------------------------------------------------------------- -# Set up the ESPResSo system and pyMBE -# ---------------------------------------------------------------------- -espresso_system = espressomd.System(box_l=[20] * 3) -pmb = pyMBE.pymbe_library(seed=42) - -# ---------------------------------------------------------------------- -# Particle definitions -# ---------------------------------------------------------------------- -pmb.define_particle(name='A', - z=0, - sigma=0.4 * pmb.units.nm, - epsilon=1 * pmb.units('reduced_energy')) - -pmb.define_particle(name='B', - z=0, - sigma=0.4 * pmb.units.nm, - epsilon=1 * pmb.units('reduced_energy')) - -pmb.define_particle(name='C', - z=0, - sigma=0.4 * pmb.units.nm, - epsilon=1 * pmb.units('reduced_energy')) - -pmb.define_particle(name='D', - z=0, - sigma=0.4 * pmb.units.nm, - epsilon=1 * pmb.units('reduced_energy')) - -# ---------------------------------------------------------------------- -# Default bond used everywhere -# ---------------------------------------------------------------------- -HARMONIC_bond_parameters = { - 'r_0': 0.4 * pmb.units.nm, - 'k': 40 * pmb.units('reduced_energy / reduced_length**2'), -} -pmb.define_default_bond(bond_type='harmonic', - bond_parameters=HARMONIC_bond_parameters) - -# ---------------------------------------------------------------------- -# Default angle used by both demos -# ---------------------------------------------------------------------- -default_angle_params = { - 'k': 50 * pmb.units('reduced_energy'), - 'phi_0': np.pi * pmb.units(''), -} -pmb.define_default_angle(angle_type="harmonic", - angle_parameters=default_angle_params) - - -def show_database(label): - print(f"\n############ {label} ############") - for pmb_type in ["particle", "residue", "molecule", "bond", "angle"]: - print(f"\n=== {pmb_type} templates ===") - print(pmb.get_templates_df(pmb_type=pmb_type)) - print(f"\n=== {pmb_type} instances ===") - print(pmb.get_instances_df(pmb_type=pmb_type)) - - -# ====================================================================== -# Demo 1: create_residue with gen_angle=True -# ====================================================================== -# Topology: central A bonded to side particles A, B, C. -# Auto-generated angles (central A has 3 neighbors): A-A-B, A-A-C, B-A-C -pmb.define_residue(name="Res_1", - central_bead="A", - side_chains=["A", "B", "C"]) - -residue_id = pmb.create_residue(name="Res_1", - espresso_system=espresso_system, - use_default_bond=True, - gen_angle=True) - -show_database("Demo 1: single residue with gen_angle=True") - -# ---------------------------------------------------------------------- -# Clean up: remove the residue (and its particles, bonds, angles) from -# both ESPResSo and the pyMBE database before running Demo 2. -# ---------------------------------------------------------------------- -pmb.delete_instances_in_system(instance_id=residue_id, - pmb_type="residue", - espresso_system=espresso_system) - -# ====================================================================== -# Demo 2: create_molecule with gen_angle=True -# ====================================================================== -# Topology: a chain of 4 Res_2 residues, each with central A and one -# side chain B. After auto-generation, angles include both intra-residue -# (A-A-B) and inter-residue (A-A-A) triplets. -pmb.define_residue(name="Res_2", - central_bead="A", - side_chains=["B"]) - -pmb.define_molecule(name="Mol_1", - residue_list=["Res_2"] * 4) - -molecule_ids = pmb.create_molecule(name="Mol_1", - number_of_molecules=1, - espresso_system=espresso_system, - backbone_vector=[1.0, 0.0, 0.0], - use_default_bond=True, - gen_angle=True) - -show_database("Demo 2: linear molecule with gen_angle=True") - -pmb.delete_instances_in_system(instance_id=molecule_ids[0], - pmb_type="molecule", - espresso_system=espresso_system) - -# ====================================================================== -# Demo 3: create_residue with a nested residue as side chain -# ====================================================================== -# Topology: D₁ where SubRes_3 = B bonded to C and D₂ -# \ -# A -# | -# B -- C -# | -# D₂ -# -# SubRes_3 (central B, side chains C and D) is passed as a side chain of -# the outer residue NestedRes_3 (central A, direct side chain D). -# D₁ and D₂ are distinct instances of particle type D. -# -# All bonded triplets and their canonical angle names: -# centered at A: D₁ -- A -- B → "B-A-D" (outer {B,D} sorted) -# centered at B: A -- B -- C → "A-B-C" -# centered at B: A -- B -- D₂ → "A-B-D" -# centered at B: C -- B -- D₂ → "C-B-D" -pmb.define_residue(name="SubRes_3", - central_bead="B", - side_chains=["C", "D"]) - -pmb.define_residue(name="NestedRes_3", - central_bead="A", - side_chains=["D", "SubRes_3"]) - -nested_residue_id = pmb.create_residue(name="NestedRes_3", - espresso_system=espresso_system, - use_default_bond=True, - gen_angle=True) - -show_database("Demo 3: nested residue with gen_angle=True") - -angle_instances_df = pmb.get_instances_df(pmb_type="angle") -particle_instances_df = pmb.get_instances_df(pmb_type="particle") -print(f"\nAngle instance count: {len(angle_instances_df)}") - -id_to_type = dict(zip(particle_instances_df["particle_id"], - particle_instances_df["name"])) -actual_triplets = set() -for _, row in angle_instances_df.iterrows(): - p1 = id_to_type[row["particle_id1"]] - center = id_to_type[row["particle_id2"]] - p3 = id_to_type[row["particle_id3"]] - outer = sorted([p1, p3]) - actual_triplets.add(f"{outer[0]}-{center}-{outer[1]}") - -print("Resolved angle triplets:", sorted(actual_triplets)) -expected_triplets = {"A-B-C", "A-B-D", "B-A-D", "C-B-D"} -assert actual_triplets == expected_triplets, ( - f"Demo 3: expected {sorted(expected_triplets)}, got {sorted(actual_triplets)}" -) -print("Demo 3: all expected angle triplets generated OK") From cac38475e594ac5f5bd33af8e8f75f1c21c9aa04 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sun, 3 May 2026 19:12:17 +0200 Subject: [PATCH 11/13] modification in setup_angular_potential.py --- samples/setup_angular_potential.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/samples/setup_angular_potential.py b/samples/setup_angular_potential.py index 80183fab..4f6bb3da 100644 --- a/samples/setup_angular_potential.py +++ b/samples/setup_angular_potential.py @@ -21,12 +21,13 @@ # This script shows three usages of `gen_angle=True`: # 1. create_residue with gen_angle=True # - builds a single residue with auto-generated angles -# 2. create_molecule with gen_angle=True -# - builds a chain of residues; angles are auto-generated across the -# full molecule, including angles that span residue boundaries -# 3. create_residue with a nested residue as side chain +# 2. create_residue with a nested residue as side chain # - builds a residue whose side chain is itself a residue; verifies # that all bonded triplets (intra- and cross-level) are generated +# 3. create_molecule with gen_angle=True +# - builds a chain of residues; angles are auto-generated across the +# full molecule, including angles that span residue boundaries + import pyMBE import numpy as np From 7d4c41d70129f918fff7deff9410de333c5d3923 Mon Sep 17 00:00:00 2001 From: Somesh Kurahatti Date: Sun, 3 May 2026 21:31:41 +0200 Subject: [PATCH 12/13] simplified build_hydrogel.py --- samples/build_hydrogel.py | 244 +++++++++++++------------------------- 1 file changed, 80 insertions(+), 164 deletions(-) diff --git a/samples/build_hydrogel.py b/samples/build_hydrogel.py index 56f430b5..df6c61db 100644 --- a/samples/build_hydrogel.py +++ b/samples/build_hydrogel.py @@ -16,10 +16,12 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . # -# Demonstrates: -# 1. Original hydrogel build (mpc=40, two bead types, no angles) with a 3D plot. -# 2. Angle-potential hydrogel (mpc=4, three bead types). -# Run with --include to also add crosslinker-adjacent angle templates. +# Builds a diamond-lattice hydrogel (bead types A, C, D) with angular potentials. +# +# By default only the chain-internal angles A-A-A and A-A-C is defined. +# Use --include_crosslinker_angles to also add A-C-D and C-D-C templates. +# Use --visualize to display a 3D plot of the lattice. +# Use --mpc to set the number of monomers per chain (default: 5). import argparse import numpy as np @@ -30,99 +32,25 @@ from pyMBE.lib.lattice import DiamondLattice -def run_simple_hydrogel(espresso_system): - """ - Build and visualize a simple hydrogel with no angular potential (mpc=40). - - Defines two bead types (C, M) and one node type, builds the diamond-lattice - hydrogel, and displays a 3D plot of the lattice. - - Args: - espresso_system (espressomd.system.System): ESPResSo system object in - which the hydrogel particles and bonds will be created. - """ - pmb = pyMBE.pymbe_library(seed=42) - mpc = 40 - NodeType = "node_type" - pmb.define_particle(name=NodeType, - sigma=0.355*pmb.units.nm, - epsilon=1*pmb.units('reduced_energy')) - BeadType1 = "C" - pmb.define_particle(name=BeadType1, - sigma=0.355*pmb.units.nm, - epsilon=1*pmb.units('reduced_energy')) - BeadType2 = "M" - pmb.define_particle(name=BeadType2, - sigma=0.355*pmb.units.nm, - epsilon=1*pmb.units('reduced_energy')) - - Res1 = "res_1" - pmb.define_residue(name=Res1, central_bead=BeadType1, side_chains=[]) - Res2 = "res_2" - pmb.define_residue(name=Res2, central_bead=BeadType2, side_chains=[]) - - residue_list = [Res1]*(mpc//2) + [Res2]*(mpc//2) - pmb.define_molecule(name="hydrogel_chain", residue_list=residue_list) - - generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2') - generic_bond_l = 0.355*pmb.units.nm - HARMONIC_parameters = {'r_0': generic_bond_l, 'k': generic_harmonic_constant} - pmb.define_bond(bond_type='harmonic', - bond_parameters=HARMONIC_parameters, - particle_pairs=[[BeadType1, BeadType1], - [BeadType1, BeadType2], - [BeadType2, BeadType2], - [NodeType, BeadType1], - [NodeType, BeadType2]]) - - diamond_lattice = DiamondLattice(mpc, generic_bond_l) - espresso_system.box_l = [diamond_lattice.box_l] * 3 - lattice_builder = pmb.initialize_lattice_builder(diamond_lattice) - - indices = diamond_lattice.indices - node_topology = [] - for index in range(len(indices)): - node_topology.append({"particle_name": NodeType, - "lattice_index": indices[index]}) - - connectivity = diamond_lattice.connectivity - node_labels = lattice_builder.node_labels - reverse_node_labels = {v: k for k, v in node_labels.items()} - connectivity_with_labels = {(reverse_node_labels[i], reverse_node_labels[j]) - for i, j in connectivity} - chain_topology = [] - for node_s, node_e in connectivity_with_labels: - chain_topology.append({'node_start': node_s, - 'node_end': node_e, - 'molecule_name': "hydrogel_chain"}) - - lattice_builder.chains = chain_topology - pmb.define_hydrogel("my_hydrogel", node_topology, chain_topology) - pmb.create_hydrogel("my_hydrogel", espresso_system) - - fig = plt.figure() - ax = fig.add_subplot(111, projection="3d") - lattice_builder.draw_lattice(ax=ax, pmb=pmb) - lattice_builder.draw_simulation_box(ax) - plt.legend(fontsize=12) - plt.show() - - -def define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_angles): +def define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_angles, mpc=5): """ - Build a hydrogel with angular potentials on the diamond lattice (mpc=4). + Build a hydrogel with angular potentials on the diamond lattice. Defines three bead types: crosslinker nodes (D), chain-end beads (C), and - internal chain beads (A). Always defines the chain-internal angle A-A-C. - When include_crosslinker_angles is True, also defines D-C-A and C-D-C angles. + internal chain beads (A). Always defines the chain-internal A-A-A and A-A-C + angular potential templates. When include_crosslinker_angles is True, also + defines the D-C-A (canonical: A-C-D) and C-D-C templates. Args: espresso_system (espressomd.system.System): ESPResSo system object in which the hydrogel particles, bonds, and angles will be created. - Its box_l is updated to match the diamond lattice dimensions. - include_crosslinker_angles (bool): If True, also define angle templates - for the crosslinker-adjacent triplets D-C-A (canonical: A-C-D) and - C-D-C in addition to the chain-internal A-A-C angle. + include_crosslinker_angles (bool): If True, also define angular potential + templates for the crosslinker-adjacent triplets D-C-A (canonical: + A-C-D) and C-D-C in addition to the chain-internal A-A-A and A-A-C + templates. + mpc (int): Number of monomers per chain. Must be >= 5 (two terminal C + beads plus at least three internal A beads to generate A-A-A angles). + Defaults to 5. Returns: tuple: @@ -130,6 +58,8 @@ def define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_ and instances for this hydrogel. - hydrogel_id (int): Assembly ID of the created hydrogel, as returned by pmb.create_hydrogel. + - lattice_builder (pyMBE.lib.lattice.LatticeBuilder): Lattice + builder instance, usable for drawing the lattice. """ pmb = pyMBE.pymbe_library(seed=42) @@ -141,7 +71,6 @@ def define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_ molecule_name = "hydrogel_chain" bond_length = 0.355 * pmb.units.nm - mpc = 4 pmb.define_particle(name=node_type, sigma=bond_length, @@ -160,12 +89,11 @@ def define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_ central_bead=internal_bead_type, side_chains=[]) pmb.define_molecule(name=molecule_name, - residue_list=[ - terminal_residue_name, - internal_residue_name, - internal_residue_name, - terminal_residue_name, - ]) + residue_list=( + [terminal_residue_name] + + [internal_residue_name] * (mpc - 2) + + [terminal_residue_name] + )) harmonic_bond_parameters = { "r_0": bond_length, @@ -184,20 +112,23 @@ def define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_ "phi_0": np.pi * pmb.units(""), } pmb.define_angular_potential(angle_type="harmonic", - angle_parameters=angle_parameters, - particle_triplets=[(chain_end_bead_type, internal_bead_type, internal_bead_type)]) + angle_parameters=angle_parameters, + particle_triplets=[ + (chain_end_bead_type, internal_bead_type, internal_bead_type), # A-A-C + (internal_bead_type, internal_bead_type, internal_bead_type), # A-A-A + ]) if include_crosslinker_angles: pmb.define_angular_potential(angle_type="harmonic", - angle_parameters=angle_parameters, - particle_triplets=[ - (node_type, chain_end_bead_type, internal_bead_type), # canonical: A-C-D - (chain_end_bead_type, node_type, chain_end_bead_type), # canonical: C-D-C - ]) + angle_parameters=angle_parameters, + particle_triplets=[ + (node_type, chain_end_bead_type, internal_bead_type), # canonical: A-C-D + (chain_end_bead_type, node_type, chain_end_bead_type), # canonical: C-D-C + ]) diamond_lattice = DiamondLattice(mpc, bond_length) espresso_system.box_l = [diamond_lattice.box_l] * 3 - pmb.initialize_lattice_builder(diamond_lattice) + lattice_builder = pmb.initialize_lattice_builder(diamond_lattice) node_topology = [ {"particle_name": node_type, "lattice_index": index} @@ -213,13 +144,14 @@ def define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_ "molecule_name": molecule_name, }) + lattice_builder.chains = chain_topology pmb.define_hydrogel(name="simple_hydrogel", node_map=node_topology, chain_map=chain_topology) hydrogel_id = pmb.create_hydrogel(name="simple_hydrogel", espresso_system=espresso_system, gen_angle=True) - return pmb, hydrogel_id + return pmb, hydrogel_id, lattice_builder def summarize_angles(pmb): @@ -244,72 +176,56 @@ def summarize_angles(pmb): return angle_templates_df, angle_instances_df, angle_counts -def run_hydrogel_with_angular_potential(label, include_crosslinker_angles, expected_angle_names, espresso_system): - """ - Clear the ESPResSo system, build one angle-hydrogel case, and assert the - expected set of canonical angle names is generated. - - Args: - label (str): Human-readable name for this case, printed as a section - header and in assertion error messages. - include_crosslinker_angles (bool): Passed directly to - define_hydrogel_with_angular_potential; controls whether - crosslinker-adjacent angle templates are included. - expected_angle_names (set of str): Set of canonical angle names - (e.g. ``{"A-A-C", "A-C-D"}``}) that must be present in the angle - instance DataFrame. An AssertionError is raised if the actual set - differs. - espresso_system (espressomd.system.System): ESPResSo system object to - clear and reuse for this case. - """ - espresso_system.part.clear() - espresso_system.bonded_inter.clear() - espresso_system.non_bonded_inter.reset() - - print(f"\n############ {label} ############") - pmb, hydrogel_id = define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_angles) - angle_templates_df, angle_instances_df, angle_counts = summarize_angles(pmb) - print(f"Hydrogel assembly id: {hydrogel_id}") - print("\nAngle templates:") - print(angle_templates_df) - print("\nAngle instance dataframe:") - print(angle_instances_df) - print("Angle instance counts:") - print(angle_counts) - - angle_names = set(angle_counts.index) - if angle_names != expected_angle_names: - raise AssertionError( - f"{label}: expected angle names {sorted(expected_angle_names)}, " - f"got {sorted(angle_names)}" - ) - print(f"{label}: OK") - - if __name__ == "__main__": parser = argparse.ArgumentParser( - description="Build hydrogels: original (mpc=40) then angle-potential (mpc=4)." + description="Build a diamond-lattice hydrogel (A/C/D beads) with angular potentials." + ) + parser.add_argument( + "--mpc", + type=int, + default=5, + help="Number of monomers per chain (default: 5).", ) parser.add_argument( - "--include", + "--include_crosslinker_angles", action="store_true", - help="Include crosslinker-adjacent angles (D-C-A and C-D-C) in the angle demo.", + help="Also define crosslinker-adjacent angle templates (A-C-D and C-D-C).", + ) + parser.add_argument( + "--visualize", + action="store_true", + help="Display a 3D plot of the hydrogel lattice after creation.", ) args = parser.parse_args() espresso_system = espressomd.System(box_l=[1.0] * 3) - # Part 1: simple hydrogel (mpc=40, no angles) - run_simple_hydrogel(espresso_system) + pmb, hydrogel_id, lattice_builder = define_hydrogel_with_angular_potential( + espresso_system, args.include_crosslinker_angles, mpc=args.mpc + ) + + if args.visualize: + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + lattice_builder.draw_lattice(ax=ax, pmb=pmb) + lattice_builder.draw_simulation_box(ax) + plt.legend(fontsize=12) + plt.show() + + angle_templates_df, angle_instances_df, angle_counts = summarize_angles(pmb) + print(f"\nHydrogel assembly id: {hydrogel_id}") + print("\nAngle templates:") + print(angle_templates_df) + print("\nAngle instance dataframe:") + print(angle_instances_df) + print("\nAngle instance counts:") + print(angle_counts) - # Part 2: angle-potential hydrogel (mpc=4); clearing is done inside run_hydrogel_with_angular_potential - if args.include: - run_hydrogel_with_angular_potential(label="Chain + crosslinker angles", - include_crosslinker_angles=True, - expected_angle_names={"A-A-C", "A-C-D", "C-D-C"}, - espresso_system=espresso_system) - else: - run_hydrogel_with_angular_potential(label="Chain angles only", - include_crosslinker_angles=False, - expected_angle_names={"A-A-C"}, - espresso_system=espresso_system) + expected_angle_names = {"A-A-A", "A-A-C", "A-C-D", "C-D-C"} if args.include_crosslinker_angles else {"A-A-A", "A-A-C"} + actual_angle_names = set(angle_counts.index) + if actual_angle_names != expected_angle_names: + raise AssertionError( + f"Expected angle names {sorted(expected_angle_names)}, " + f"got {sorted(actual_angle_names)}" + ) + print("\nOK: all expected angular potential triplets generated.") From bd2448a53c25d7305c9fb05648f29548c10d66d7 Mon Sep 17 00:00:00 2001 From: Pablo Date: Sun, 3 May 2026 21:43:30 +0200 Subject: [PATCH 13/13] update changelog --- CHANGELOG.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a1e3fad..5b513a04 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ## Added +- Support to setup angular potentials with pyMBE. All flexible pyMBE templates now support angula potentials: hydrogels, molecules, peptides and residues (including residues with nested residues). (#150) +- Sample scripts and tests for the new functionality. (#150) +- Template and instance `Angle` to store information about angular potentials in the pyMBE database. (#150) +- Functions `define_angle` and ``define_default_angle` to define a templates of angular potentials in pyMBE. (#150) +- Function `create_angle` to create instances of angular potential in pyMBE. (#150) +- Function `get_angle_template` to retrieve a template of an angle potential in the pyMBE database. (#150) - Introduced a canonical pyMBE database backend replacing the previous monolithic Pandas DataFrame storage approach. This lays the foundation for more robust, extensible, and normalized data handling across pyMBE. (#147) - Added support to define reaction templates in the pyMBE database. (#147) - Utility functions to cast information about templates and instances in the pyMBE database into pandas dataframe `pmb.get_templates_df`, `pmb.get_instances_df` and `pmb.get_reactions_df`. (#147)