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) diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index f2bd213b..22746bc0 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 @@ -225,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. @@ -241,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. @@ -260,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: @@ -279,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"], @@ -296,25 +301,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_angular_potential(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. @@ -901,7 +964,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' @@ -915,6 +978,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. """ @@ -925,6 +993,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 @@ -942,21 +1011,65 @@ 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)) 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 +1208,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 +1364,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 +1488,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 +1574,246 @@ def define_default_bond(self, bond_type, bond_parameters): bond_type=bond_type) tpl.name = "default" self.db._register_template(tpl) - + + def define_angular_potential(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_angular_potential(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_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. + + 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) + + 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 + + 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] + try: + self.create_angular_potential(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/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 f4b3f404..f8ad826a 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 @@ -88,8 +90,9 @@ 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= {} 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_hydrogel.py b/samples/build_hydrogel.py index b5ab00a5..df6c61db 100644 --- a/samples/build_hydrogel.py +++ b/samples/build_hydrogel.py @@ -15,94 +15,217 @@ # # You should have received a copy of the GNU General Public License # along with this program. If not, see . +# +# 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 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 define_hydrogel_with_angular_potential(espresso_system, include_crosslinker_angles, mpc=5): + """ + 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 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. + 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: + - 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. + - lattice_builder (pyMBE.lib.lattice.LatticeBuilder): Lattice + builder instance, usable for drawing the lattice. + """ + 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 + + 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] * (mpc - 2) + + [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_angular_potential(angle_type="harmonic", + 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 + ]) + + diamond_lattice = DiamondLattice(mpc, bond_length) + espresso_system.box_l = [diamond_lattice.box_l] * 3 + lattice_builder = 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, + }) + + 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, lattice_builder + + +def summarize_angles(pmb): + """ + 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 + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + 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_crosslinker_angles", + action="store_true", + 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) + + 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) + + 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.") diff --git a/samples/setup_angular_potential.py b/samples/setup_angular_potential.py new file mode 100644 index 00000000..4f6bb3da --- /dev/null +++ b/samples/setup_angular_potential.py @@ -0,0 +1,201 @@ +# +# 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_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 +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/CTestTestfile.cmake b/testsuite/CTestTestfile.cmake index aa66fb39..63caeaea 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) @@ -75,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/angle_tests.py b/testsuite/angle_tests.py new file mode 100644 index 00000000..14f97f8f --- /dev/null +++ b/testsuite/angle_tests.py @@ -0,0 +1,790 @@ +# +# Copyright (C) 2026 pyMBE-dev team +# +# This file is part of pyMBE. +# +# pyMBE is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# pyMBE is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +# Import 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_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'), + } + + 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_angular_potential(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_angular_potential(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_angular_potential(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_angular_potential(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_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_angular_potential(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_angular_potential(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_angular_potential(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. + """ + 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_angular_potential(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_angular_potential(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_angular_potential(angle_type="cosine", + angle_parameters=self.cosine_angle_params) + + # Define specific angle for A-B-C with harmonic parameters + pmb.define_angular_potential(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_angular_potential(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_angular_potential, pmb.define_default_angular_potential]: + 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_angular_potential: + 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_angular_potential: + 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_angular_potential: + input_parameters["particle_triplets"] = [('A', 'B', 'C')] + + np.testing.assert_raises(ValueError, callback, **input_parameters) + + # 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_angular_potential, + **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_angular_potential(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_angular_potential(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_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) + + # 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_angular_potential(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_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) + + # 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_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_angular_potential(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_angular_potential(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_angular_potential(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() 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/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..632c1f74 --- /dev/null +++ b/testsuite/hydrogel_builder_with_angles.py @@ -0,0 +1,229 @@ +# +# 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 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. + + 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" + 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_angular_potential(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_angular_potential(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 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 + for angle in pmb_local.db.get_instances(pmb_type="angle").values() + ) + + +def expected_crosslinker_angle_counts(diamond_lattice_local): + """ + 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 + 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() - diff --git a/testsuite/test_io_database.py b/testsuite/test_io_database.py index 7c03c522..186eaa5c 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_angular_potential(angle_type="harmonic", + angle_parameters=parameters1, + particle_triplets=[("Z", "X", "Z"), + ("X", "Z", "Y")]) + pmb.define_default_angular_potential(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 @@ -437,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, @@ -461,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_angular_potential(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" @@ -497,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" @@ -539,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"), @@ -580,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"), @@ -755,4 +795,3 @@ def test_load_database_missing_folder(self): if __name__ == '__main__': ut.main() -