From ddc05b5846acc6c6b3f246f9e922bc88fd4a7f96 Mon Sep 17 00:00:00 2001 From: Raj Kumar Rajaram Baskaran <12559850+rajkumarrb@users.noreply.github.com> Date: Wed, 2 Sep 2026 14:07:22 +0100 Subject: [PATCH 1/4] Add support for dihedral potentials Completes the last unchecked task of pyMBE-dev/pyMBE#143 (angles were merged in #150; dihedrals are the remaining bonded-potential scope). The motivating case is Martini 2.1 peptide backbone dihedrals. New: - DihedralTemplate / DihedralInstance in pyMBE/storage. The template holds a list of terms, so GROMACS-style stacked multiplicities on the same quadruplet register as one Dihedral object per term. - pmb.define_dihedral_potential / define_default_dihedral_potential / create_dihedral_potential / get_dihedral_template. - _generate_dihedrals_for_entity and a gen_dihedral kwarg on create_residue / create_molecule / create_hydrogel, mirroring gen_angle. - espresso_engine _add_dihedral / _get_dihedral_instance / _create_dihedral_instance, wired into add_instances_to_engine. - storage manager + io: register / update / delete / templates-df / save / load branches for the dihedral type. - testsuite/dihedral_tests.py, samples/setup_dihedral_potential.py. The canonical template key is invariant under full reversal of the quadruplet. The "periodic" type maps onto the ESPResSo Dihedral bond V(phi) = k*[1 - cos(n*phi - phi_0)]; "tabulated" is reserved but raises NotImplementedError. GROMACS proper dihedral type 1, V(phi) = k*[1 + cos(n*phi - phi_s)], is recovered with phi_0 = phi_s + pi (documented in the docstrings and the sample, not applied implicitly). --- CHANGELOG.md | 7 + pyMBE/pyMBE.py | 269 +++++++++- pyMBE/simulation_builder/espresso_engine.py | 97 +++- pyMBE/storage/instances/dihedral.py | 66 +++ pyMBE/storage/io.py | 57 ++- pyMBE/storage/manager.py | 48 +- pyMBE/storage/templates/dihedral.py | 142 ++++++ samples/setup_dihedral_potential.py | 180 +++++++ testsuite/CTestTestfile.cmake | 1 + testsuite/database_unit_tests.py | 37 +- testsuite/dihedral_tests.py | 523 ++++++++++++++++++++ testsuite/test_io_database.py | 66 +++ 12 files changed, 1478 insertions(+), 15 deletions(-) create mode 100644 pyMBE/storage/instances/dihedral.py create mode 100644 pyMBE/storage/templates/dihedral.py create mode 100644 samples/setup_dihedral_potential.py create mode 100644 testsuite/dihedral_tests.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a15004b..019daec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ## Added +- Support to setup dihedral (torsion) potentials with pyMBE, completing the bonded-potential support tracked in #143. All flexible pyMBE templates support dihedral potentials: molecules, peptides, residues (including residues with nested residues) and the internal chains of hydrogels. (#156) +- Template and instance `Dihedral` to store information about dihedral potentials in the pyMBE database, including support for multi-term (stacked multiplicity) dihedrals on the same quadruplet. (#156) +- Functions `define_dihedral_potential` and `define_default_dihedral_potential` to define templates of dihedral potentials in pyMBE. (#156) +- Function `create_dihedral_potential` to create instances of dihedral potentials in pyMBE. (#156) +- Function `get_dihedral_template` to retrieve a template of a dihedral potential in the pyMBE database. (#156) +- `gen_dihedral` keyword argument in `create_residue`, `create_molecule` and `create_hydrogel` to auto-generate dihedral potentials from the bond topology. (#156) +- Sample script `samples/setup_dihedral_potential.py` and test suite `testsuite/dihedral_tests.py` for the new functionality. (#156) - Simulation-engine abstraction layer to decouple pyMBE molecular-building logic from simulation-engine-specific implementation details. (#151) - Added `EspressoEngine` to delegate ESPResSo-specific operations, including reaction methods, setup utilities, and particle/bond/angle insertion. (#151) - Added `add_instances_to_engine()` to transfer pyMBE particles, bonds, and angles to the selected simulation engine. (#151) diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index 5b0f8ca..e1d582c 100644 --- a/pyMBE/pyMBE.py +++ b/pyMBE/pyMBE.py @@ -39,6 +39,7 @@ 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.dihedral import DihedralTemplate ## Instances from pyMBE.storage.instances.particle import ParticleInstance from pyMBE.storage.instances.residue import ResidueInstance @@ -47,6 +48,7 @@ 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.dihedral import DihedralInstance from pyMBE.storage.instances.hydrogel import HydrogelInstance from pyMBE.simulation_builder.espresso_engine import EspressoSimulation @@ -193,7 +195,7 @@ def _check_pka_set(self, pka_set): if required_key not in pka_entry: raise ValueError(f'missing a required key "{required_key}" in entry "{pka_name}" of pka_set ("{pka_entry}")') - def _create_hydrogel_chain(self, hydrogel_chain, nodes,box_l, use_default_bond=False, gen_angle=False): + def _create_hydrogel_chain(self, hydrogel_chain, nodes,box_l, use_default_bond=False, gen_angle=False, gen_dihedral=False): """ Creates a chain between two nodes of a hydrogel. @@ -211,6 +213,10 @@ def _create_hydrogel_chain(self, hydrogel_chain, nodes,box_l, use_default_bond=F chain molecule. Junction angles near the hydrogel crosslinkers are handled separately at the hydrogel level. + gen_dihedral ('bool', optional): + If True, generate the dihedral potentials internal to the created + chain molecule. + Return: ('int'): molecule_id of the created hydrogel chian. @@ -265,7 +271,8 @@ def _create_hydrogel_chain(self, hydrogel_chain, nodes,box_l, use_default_bond=F backbone_vector=np.array(backbone_vector)/l0, use_default_bond=use_default_bond, reverse_residue_order=reverse_residue_order, - gen_angle=gen_angle)[0] + gen_angle=gen_angle, + gen_dihedral=gen_dihedral)[0] chain_pids = self.db._find_instance_ids_by_attribute(pmb_type="particle", attribute="molecule_id", value=mol_id) @@ -833,7 +840,7 @@ def create_counterions(self, object_name, cation_name, anion_name, box_l): return counterion_number - def create_hydrogel(self, name, box_l, use_default_bond=False, gen_angle=False): + def create_hydrogel(self, name, box_l, use_default_bond=False, gen_angle=False, gen_dihedral=False): """ Creates a hydrogel in espresso_system using a pyMBE hydrogel template given by 'name' @@ -851,6 +858,11 @@ def create_hydrogel(self, name, box_l, use_default_bond=False, gen_angle=False): chains and, when explicitly defined, for all crosslinker-adjacent triplets. Defaults to False. + gen_dihedral ('bool', optional): + If True, generate dihedral potentials for the internal hydrogel + chains. Crosslinker-adjacent dihedrals are not generated. + Defaults to False. + Returns: ('int'): id of the hydrogel instance created. """ @@ -881,6 +893,7 @@ def create_hydrogel(self, name, box_l, use_default_bond=False, gen_angle=False): box_l=box_l, use_default_bond=use_default_bond, gen_angle=gen_angle, + gen_dihedral=gen_dihedral, ) self.db._update_instance(instance_id=molecule_id, pmb_type="molecule", @@ -939,7 +952,7 @@ def create_hydrogel(self, name, box_l, use_default_bond=False, gen_angle=False): return assembly_id - def create_molecule(self, name, number_of_molecules, box_l, list_of_first_residue_positions=None, backbone_vector=None, use_default_bond=False, reverse_residue_order = False, gen_angle=False): + def create_molecule(self, name, number_of_molecules, box_l, list_of_first_residue_positions=None, backbone_vector=None, use_default_bond=False, reverse_residue_order = False, gen_angle=False, gen_dihedral=False): """ Creates instances of a given molecule template name into ESPResSo. @@ -1081,6 +1094,10 @@ def create_molecule(self, name, number_of_molecules, box_l, list_of_first_residu self._generate_angles_for_entity( entity_id=molecule_id, entity_id_col='molecule_id') + if gen_dihedral: + self._generate_dihedrals_for_entity( + entity_id=molecule_id, + entity_id_col='molecule_id') first_residue = True pos_index+=1 molecule_ids.append(molecule_id) @@ -1230,7 +1247,7 @@ def create_protein(self, name, number_of_proteins, box_l, topology_dict): mol_ids.append(molecule_id) return mol_ids - def create_residue(self, name, box_l, central_bead_position=None,use_default_bond=False, backbone_vector=None, gen_angle=False): + def create_residue(self, name, box_l, central_bead_position=None,use_default_bond=False, backbone_vector=None, gen_angle=False, gen_dihedral=False): """ Creates a residue into ESPResSo. @@ -1361,6 +1378,10 @@ def create_residue(self, name, box_l, central_bead_position=None,use_default_bon self._generate_angles_for_entity( entity_id=residue_id, entity_id_col="residue_id") + if gen_dihedral: + self._generate_dihedrals_for_entity( + entity_id=residue_id, + entity_id_col="residue_id") return residue_id def define_bond(self, bond_type, bond_parameters, particle_pairs): @@ -1625,6 +1646,244 @@ def _generate_angles_for_entity(self, entity_id, entity_id_col): # No angle template defined for this triplet — skip continue + def _build_dihedral_terms(self, dihedral_type, dihedral_parameters): + """ + Validate user input for a dihedral potential and cast it into the list + of ``PintQuantity`` terms stored in a `DihedralTemplate`. + + Args: + dihedral_type ('str'): + Type of dihedral potential. Supported: "periodic". + + dihedral_parameters ('dict' or 'list[dict]'): + Parameters of the dihedral potential. GROMACS-style force fields + stack several terms (with different multiplicities) on the same + quadruplet, so a list of term dictionaries is accepted. A single + dictionary is treated as a one-term dihedral. For the "periodic" + type each term must contain: + - "k" ('pint.Quantity'): Torsion stiffness with dimensions of energy. + - "phi_0" ('pint.Quantity'): Equilibrium phase, dimensionless (radians). + - "mult" ('int'): Multiplicity (number of minima). + + Returns: + ('list[dict[str, PintQuantity]]'): One dictionary per dihedral term. + """ + valid_dihedral_types = ["periodic", "tabulated"] + if dihedral_type not in valid_dihedral_types: + raise NotImplementedError(f"Dihedral potential type '{dihedral_type}' currently not implemented in pyMBE, accepted types are {valid_dihedral_types}") + if dihedral_type == "tabulated": + raise NotImplementedError("Tabulated dihedral potentials are not supported yet") + + if isinstance(dihedral_parameters, dict): + term_list = [dihedral_parameters] + else: + term_list = list(dihedral_parameters) + if not term_list: + raise ValueError("At least one dihedral term must be provided") + + terms = [] + for term in term_list: + if "k" not in term: + raise ValueError("Magnitude of the dihedral potential (k) is missing") + if "phi_0" not in term: + raise ValueError("Equilibrium phase of the dihedral potential (phi_0) is missing") + if "mult" not in term: + raise ValueError("Multiplicity of the dihedral potential (mult) is missing") + mult = term["mult"] + if not isinstance(mult, self.units.Quantity): + mult = mult * self.units("") + terms.append({"k": PintQuantity.from_quantity(q=term["k"], + expected_dimension="energy", + ureg=self.units), + "phi_0": PintQuantity.from_quantity(q=term["phi_0"], + expected_dimension="dimensionless", + ureg=self.units), + "mult": PintQuantity.from_quantity(q=mult, + expected_dimension="dimensionless", + ureg=self.units)}) + return terms + + def define_dihedral_potential(self, dihedral_type, dihedral_parameters, particle_quadruplets): + """ + Defines dihedral potential templates for each particle quadruplet in `particle_quadruplets`. + + Args: + dihedral_type ('str'): + Type of dihedral potential. Supported: "periodic", which maps onto + the ESPResSo ``Dihedral`` bond ``V(phi) = k*[1 - cos(n*phi - phi_0)]``. + + dihedral_parameters ('dict' or 'list[dict]'): + Parameters of the dihedral potential, or a list of such + dictionaries to register several terms (different multiplicities) + on the same quadruplet. For the "periodic" type each term must + contain: + - "k" ('pint.Quantity'): Torsion stiffness with dimensions of energy. + - "phi_0" ('pint.Quantity'): Equilibrium phase, dimensionless (radians). + - "mult" ('int'): Multiplicity (number of minima). + + particle_quadruplets ('list[tuple[str,str,str,str]]'): + List of (particle1, particle2, particle3, particle4) name quadruplets. + + Notes: + GROMACS proper dihedral type 1, ``V(phi) = k*[1 + cos(n*phi - phi_s)]``, + is recovered by passing ``k`` unchanged, ``mult = n`` and + ``phi_0 = phi_s + pi`` (with ``phi_s`` converted from degrees to radians). + """ + terms = self._build_dihedral_terms(dihedral_type=dihedral_type, + dihedral_parameters=dihedral_parameters) + dihedral_names = [] + for name1, name2, name3, name4 in particle_quadruplets: + tpl = DihedralTemplate(particle_name1=name1, + particle_name2=name2, + particle_name3=name3, + particle_name4=name4, + terms=terms, + dihedral_type=dihedral_type) + tpl._make_name() + if tpl.name in dihedral_names: + raise RuntimeError(f"Dihedral {tpl.name} has already been defined, please check the list of particle quadruplets") + dihedral_names.append(tpl.name) + self.db._register_template(tpl) + + def define_default_dihedral_potential(self, dihedral_type, dihedral_parameters): + """ + Defines a dihedral template as a "default" template in the pyMBE database. + + Args: + dihedral_type ('str'): + Type of dihedral potential. Supported: "periodic". + + dihedral_parameters ('dict' or 'list[dict]'): + Parameters of the dihedral potential (k, phi_0, mult), or a list + of such dictionaries for a multi-term dihedral. + """ + terms = self._build_dihedral_terms(dihedral_type=dihedral_type, + dihedral_parameters=dihedral_parameters) + tpl = DihedralTemplate(terms=terms, + dihedral_type=dihedral_type) + tpl.name = "default" + self.db._register_template(tpl) + + def create_dihedral_potential(self, particle_id1, particle_id2, particle_id3, particle_id4, use_default_dihedral=False): + """ + Creates a dihedral between four particle instances in an ESPResSo system + and registers it in the pyMBE database. + + Args: + particle_id1 ('int'): ID of the first particle in the quadruplet. + particle_id2 ('int'): ID of the second particle in the quadruplet. + particle_id3 ('int'): ID of the third particle in the quadruplet. + particle_id4 ('int'): ID of the fourth particle in the quadruplet. + use_default_dihedral ('bool', optional): If True, use the default dihedral 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) + particle_inst_4 = self.db.get_instance(pmb_type="particle", instance_id=particle_id4) + + # Verify that the i-j, j-k and k-l bonds of the quadruplet exist + bond_instances = self.db.get_instances(pmb_type="bond") + bonded_pairs = set() + for bond in bond_instances.values(): + bonded_pairs.add(frozenset([bond.particle_id1, bond.particle_id2])) + chain = [(particle_id1, particle_id2), + (particle_id2, particle_id3), + (particle_id3, particle_id4)] + for pid_a, pid_b in chain: + if frozenset([pid_a, pid_b]) not in bonded_pairs: + raise ValueError(f"Cannot create dihedral: no bond exists between particle {pid_a} and particle {pid_b}.") + + dihedral_tpl = self.get_dihedral_template(name1=particle_inst_1.name, + name2=particle_inst_2.name, + name3=particle_inst_3.name, + name4=particle_inst_4.name, + use_default_dihedral=use_default_dihedral) + dihedral_id = self.db._propose_instance_id(pmb_type="dihedral") + pmb_dihedral_instance = DihedralInstance(dihedral_id=dihedral_id, + name=dihedral_tpl.name, + particle_id1=particle_id1, + particle_id2=particle_id2, + particle_id3=particle_id3, + particle_id4=particle_id4) + self.db._register_instance(instance=pmb_dihedral_instance) + + def get_dihedral_template(self, name1, name2, name3, name4, use_default_dihedral=False): + """ + Retrieves a dihedral template connecting four particle templates. + + Args: + name1 ('str'): Name of the first particle in the quadruplet. + name2 ('str'): Name of the second particle in the quadruplet. + name3 ('str'): Name of the third particle in the quadruplet. + name4 ('str'): Name of the fourth particle in the quadruplet. + use_default_dihedral ('bool', optional): If True, fall back to the default dihedral template. + + Returns: + ('DihedralTemplate'): The matching dihedral template. + """ + dihedral_key = DihedralTemplate.make_dihedral_key(name1, name2, name3, name4) + try: + return self.db.get_template(name=dihedral_key, pmb_type="dihedral") + except ValueError: + pass + + if use_default_dihedral: + return self.db.get_template(name="default", pmb_type="dihedral") + + raise ValueError(f"No dihedral template found for '{name1}-{name2}-{name3}-{name4}', and default dihedrals are deactivated.") + + def _generate_dihedrals_for_entity(self, entity_id, entity_id_col): + """ + Auto-generates dihedrals from bond topology for an entity (molecule or residue). + + For each bond ``(j, k)`` internal to the entity, this method forms the + candidate quadruplets ``(i, j, k, l)`` with ``i`` a bonded neighbour of + ``j`` and ``l`` a bonded neighbour of ``k``, and applies any matching + dihedral potential. + + Args: + entity_id ('int'): The molecule_id or residue_id to generate dihedrals 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) + + seen_quadruplets = set() + for j in particle_ids: + for k in sorted(neighbors[j]): + # Iterate over each internal bond (j, k) once + for i in sorted(neighbors[j] - {k}): + for l in sorted(neighbors[k] - {j}): + if i == l: + continue + quad = (i, j, k, l) + if quad in seen_quadruplets or (l, k, j, i) in seen_quadruplets: + continue + seen_quadruplets.add(quad) + try: + self.create_dihedral_potential(particle_id1=i, + particle_id2=j, + particle_id3=k, + particle_id4=l, + use_default_dihedral=True) + except ValueError: + # No dihedral template defined for this quadruplet — skip + continue + def define_hydrogel(self, name, node_map, chain_map): """ Defines a hydrogel template in the pyMBE database. diff --git a/pyMBE/simulation_builder/espresso_engine.py b/pyMBE/simulation_builder/espresso_engine.py index 68fe830..dbae054 100644 --- a/pyMBE/simulation_builder/espresso_engine.py +++ b/pyMBE/simulation_builder/espresso_engine.py @@ -57,7 +57,29 @@ def _add_angle(self,particle_id1,particle_id2,particle_id3, angle_inst): pmb_type='angle', attribute='added_to_engine', value=True) - + + def _add_dihedral(self,particle_id1,particle_id2,particle_id3,particle_id4, dihedral_inst): + """ + Helper function to add dihedral instances to ESPResSo. + + Args: + particle_id1 (int): pid of the particle instance 1 that composes the dihedral + particle_id2 (int): pid of the particle instance 2 that composes the dihedral + particle_id3 (int): pid of the particle instance 3 that composes the dihedral + particle_id4 (int): pid of the particle instance 4 that composes the dihedral + dihedral_inst (pmb.DihedralInstance): dataclass containing information of a dihedral instance + """ + dihedral_tpl=self.db.get_template(name=dihedral_inst.name, + pmb_type="dihedral") + espresso_dihedral_insts=self._get_dihedral_instance(dihedral_template=dihedral_tpl) + + for espresso_dihedral_inst in espresso_dihedral_insts: + self.espresso_system.part.by_id(particle_id1).add_bond((espresso_dihedral_inst, particle_id2, particle_id3, particle_id4)) + self.db._update_instance(instance_id=dihedral_inst.dihedral_id, + pmb_type='dihedral', + attribute='added_to_engine', + value=True) + def _add_bond(self,particle_id1,particle_id2,bond_inst): """ Helper function to add bond instances to ESPResSo. @@ -248,7 +270,57 @@ def _create_angle_instance(self, angle_type, angle_parameters): elif angle_type == "harmonic_cosine": return interactions.AngleCossquare(bend=k, phi0=phi_0) - + def _get_dihedral_instance(self,dihedral_template): + """ + Retrieve or create the ESPResSo dihedral interactions for a given dihedral template. + + A dihedral template can hold several terms (different multiplicities on + the same quadruplet), so this returns a list of ESPResSo interaction + objects, one per term. + + Args: + dihedral_template ('DihedralTemplate'): The dihedral template to use. + + Returns: + ('list[espressomd.interactions.BondedInteraction]'): The ESPResSo dihedral interaction objects. + """ + if dihedral_template.name in self.db.espresso_dihedral_instances: + return self.db.espresso_dihedral_instances[dihedral_template.name] + + dihedral_insts = self._create_dihedral_instance(dihedral_type=dihedral_template.dihedral_type, + dihedral_terms=dihedral_template.get_terms(self.units)) + self.db.espresso_dihedral_instances[dihedral_template.name] = dihedral_insts + for dihedral_inst in dihedral_insts: + self.espresso_system.bonded_inter.add(dihedral_inst) + + return dihedral_insts + + def _create_dihedral_instance(self, dihedral_type, dihedral_terms): + """ + Creates the ESPResSo dihedral interaction objects for a list of terms. + + Args: + dihedral_type ('str'): Type of dihedral potential ("periodic"). + dihedral_terms ('list[dict]'): List of term dictionaries (k, phi_0, mult). + + Returns: + ('list[espressomd.interactions.BondedInteraction]'): The ESPResSo dihedral interaction objects. + """ + from espressomd import interactions + + if dihedral_type == "periodic": + dihedral_insts = [] + for term in dihedral_terms: + k = term["k"].m_as("reduced_energy") + phi_0 = float(term["phi_0"].magnitude) + mult = int(round(float(term["mult"].magnitude))) + dihedral_insts.append(interactions.Dihedral(bend=k, mult=mult, phase=phi_0)) + return dihedral_insts + if dihedral_type == "tabulated": + raise NotImplementedError("tabulated dihedrals not yet supported") + raise NotImplementedError(f"Dihedral potential type '{dihedral_type}' currently not implemented in pyMBE") + + def _get_particle_ids_in_espresso(self): """ Gets a list of all the particles_id in espresso_system instance. @@ -1354,8 +1426,8 @@ def update_particle_id(self, old_pid, new_pid): def add_instances_to_engine(self): """ - This method adds the set of particles instances and bond instances and angle instances - that are not present in the pymbe data base + This method adds the set of particles instances, bond instances, angle instances + and dihedral instances that are not present in the pymbe data base """ missing_particle_ids=self.db._find_instance_ids_by_attribute(pmb_type='particle', @@ -1413,4 +1485,21 @@ def add_instances_to_engine(self): particle_id2, particle_id3, angle_instance) + + missing_dihedral_ids=self.db._find_instance_ids_by_attribute(pmb_type='dihedral', + attribute='added_to_engine', + value=False) + if len(missing_dihedral_ids)>0: + for id in missing_dihedral_ids: + dihedral_instance=self.db.get_instance(pmb_type='dihedral', + instance_id=id) + particle_id1=dihedral_instance.particle_id1 + particle_id2=dihedral_instance.particle_id2 + particle_id3=dihedral_instance.particle_id3 + particle_id4=dihedral_instance.particle_id4 + self._add_dihedral(particle_id1, + particle_id2, + particle_id3, + particle_id4, + dihedral_instance) diff --git a/pyMBE/storage/instances/dihedral.py b/pyMBE/storage/instances/dihedral.py new file mode 100644 index 0000000..aa42bad --- /dev/null +++ b/pyMBE/storage/instances/dihedral.py @@ -0,0 +1,66 @@ +# +# 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 DihedralInstance(PMBBaseModel): + """ + Instance representation of a dihedral (torsion) between four particles. + + Attributes: + pmb_type ('Literal["dihedral"]'): + Fixed identifier set to ``"dihedral"`` for all dihedral instances. + + dihedral_id ('int'): + Unique non-negative integer identifying this dihedral instance. + + name ('str'): + Name of the dihedral template from which this instance was created. + + particle_id1 ('int'): + ID of the first particle in the quadruplet. + + particle_id2 ('int'): + ID of the second particle in the quadruplet. + + particle_id3 ('int'): + ID of the third particle in the quadruplet. + + particle_id4 ('int'): + ID of the fourth particle in the quadruplet. + + added_to_engine ('bool'): + Whether the dihedral has already been transferred to the simulation engine. + """ + pmb_type: Literal["dihedral"] = "dihedral" + dihedral_id: int + name: str + particle_id1: int + particle_id2: int + particle_id3: int + particle_id4: int + added_to_engine: bool = False + + @validator("dihedral_id", "particle_id1", "particle_id2", "particle_id3", "particle_id4") + 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 0469ce1..ae55f62 100644 --- a/pyMBE/storage/io.py +++ b/pyMBE/storage/io.py @@ -20,7 +20,7 @@ import os import json from pathlib import Path -from typing import Any, Dict +from typing import Any, Dict, List import pandas as pd import numpy as np import logging @@ -31,11 +31,13 @@ from pyMBE.storage.templates.molecule import MoleculeTemplate from pyMBE.storage.templates.bond import BondTemplate from pyMBE.storage.templates.angle import AngleTemplate +from pyMBE.storage.templates.dihedral import DihedralTemplate 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.instances.dihedral import DihedralInstance from pyMBE.storage.reactions.reaction import Reaction, ReactionParticipant from pyMBE.storage.templates.peptide import PeptideTemplate from pyMBE.storage.instances.peptide import PeptideInstance @@ -126,7 +128,7 @@ def _load_database_csv(db, folder): Notes: - PintQuantity objects are reconstructed from their dictionary representation. - - Supports particle, residue, molecule, peptide, protein, bond, angle, and hydrogel types. + - Supports particle, residue, molecule, peptide, protein, bond, angle, dihedral, and hydrogel types. """ folder = Path(folder) if not folder.exists(): @@ -136,8 +138,9 @@ def _load_database_csv(db, folder): "particle_state", "residue", "molecule", - "bond", + "bond", "angle", + "dihedral", "peptide", "protein", "hydrogel", @@ -238,6 +241,24 @@ def _load_database_csv(db, folder): side_particle2=None if side_particle2 == "" else side_particle2, parameters=parameters) templates[tpl.name] = tpl + elif pmb_type == "dihedral": + terms_raw = _decode(row.get("terms", "")) or [] + terms: List[Dict[str, Any]] = [] + for term_raw in terms_raw: + term: Dict[str, Any] = {} + for k, v in term_raw.items(): + if isinstance(v, dict) and {"magnitude", "units", "dimension"}.issubset(v.keys()): + term[k] = PintQuantity.from_dict(v) + terms.append(term) + particle_names = [row.get(f"particle_name{i}", "") or "" for i in range(1, 5)] + tpl = DihedralTemplate(name=row["name"], + dihedral_type=row.get("dihedral_type", ""), + particle_name1=None if particle_names[0] == "" else particle_names[0], + particle_name2=None if particle_names[1] == "" else particle_names[1], + particle_name3=None if particle_names[2] == "" else particle_names[2], + particle_name4=None if particle_names[3] == "" else particle_names[3], + terms=terms) + 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 [] @@ -340,6 +361,15 @@ def _load_database_csv(db, folder): particle_id3=int(row["particle_id3"]), added_to_engine=row["added_to_engine"]) instances[inst.angle_id] = inst + elif pmb_type == "dihedral": + inst = DihedralInstance(name=row["name"], + dihedral_id=int(row["dihedral_id"]), + particle_id1=int(row["particle_id1"]), + particle_id2=int(row["particle_id2"]), + particle_id3=int(row["particle_id3"]), + particle_id4=int(row["particle_id4"]), + added_to_engine=row["added_to_engine"]) + instances[inst.dihedral_id] = inst elif pmb_type == "hydrogel": inst = HydrogelInstance(name=row["name"], assembly_id=int(row["assembly_id"])) @@ -444,6 +474,18 @@ def _save_database_csv(db, folder): "side_particle2": tpl.side_particle2, "angle_type": tpl.angle_type, "parameters": _encode(params_serial)}) + elif pmb_type == "dihedral" and isinstance(tpl, DihedralTemplate): + terms_serial = [] + for term in tpl.terms: + terms_serial.append({k: v.to_dict() for k, v in term.items() + if isinstance(v, PintQuantity)}) + rows.append({"name": tpl.name, + "particle_name1": tpl.particle_name1, + "particle_name2": tpl.particle_name2, + "particle_name3": tpl.particle_name3, + "particle_name4": tpl.particle_name4, + "dihedral_type": tpl.dihedral_type, + "terms": _encode(terms_serial)}) # HYDROGEL TEMPLATE elif pmb_type == "hydrogel" and isinstance(tpl, HydrogelTemplate): rows.append({"name": tpl.name, @@ -524,6 +566,15 @@ def _save_database_csv(db, folder): "particle_id2": int(inst.particle_id2), "particle_id3": int(inst.particle_id3), "added_to_engine":inst.added_to_engine,}) + elif pmb_type == "dihedral" and isinstance(inst, DihedralInstance): + rows.append({"pmb_type": pmb_type, + "name": inst.name, + "dihedral_id": int(inst.dihedral_id), + "particle_id1": int(inst.particle_id1), + "particle_id2": int(inst.particle_id2), + "particle_id3": int(inst.particle_id3), + "particle_id4": int(inst.particle_id4), + "added_to_engine":inst.added_to_engine,}) elif pmb_type == "hydrogel" and isinstance(inst, HydrogelInstance): rows.append({"pmb_type": pmb_type, "name": inst.name, diff --git a/pyMBE/storage/manager.py b/pyMBE/storage/manager.py index 1ced265..f32ba10 100644 --- a/pyMBE/storage/manager.py +++ b/pyMBE/storage/manager.py @@ -27,11 +27,13 @@ from pyMBE.storage.templates.molecule import MoleculeTemplate from pyMBE.storage.templates.bond import BondTemplate from pyMBE.storage.templates.angle import AngleTemplate +from pyMBE.storage.templates.dihedral import DihedralTemplate 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.instances.dihedral import DihedralInstance from pyMBE.storage.reactions.reaction import Reaction from pyMBE.storage.templates.peptide import PeptideTemplate from pyMBE.storage.instances.peptide import PeptideInstance @@ -92,9 +94,10 @@ def __init__(self,units): "peptide", "protein"] self._assembly_like_types = ["hydrogel"] - self._pmb_types = ["particle", "residue", "angle"] + self._molecule_like_types + self._assembly_like_types + self._pmb_types = ["particle", "residue", "angle", "dihedral"] + self._molecule_like_types + self._assembly_like_types self.espresso_bond_instances= {} self.espresso_angle_instances= {} + self.espresso_dihedral_instances= {} def _collect_particle_templates(self, name, pmb_type): """ @@ -203,6 +206,28 @@ def _delete_angles_of_particle(self, pid): if "angle" in self._instances and not self._instances["angle"]: del self._instances["angle"] + def _delete_dihedrals_of_particle(self, pid): + """ + Delete all dihedral instances involving a given particle instance. + + Args: + pid ('int'): + The particle ID whose associated dihedrals should be deleted. + + Notes: + - If no `"dihedral"` instances are present in the database, the method + exits immediately. + - This method does not raise errors if no dihedrals involve the particle. + - It is intended for internal use by cascade-deletion routines. + """ + if "dihedral" not in self._instances: + return + dihedrals_to_delete = [d_id for d_id, d in list(self._instances["dihedral"].items()) if pid in (d.particle_id1, d.particle_id2, d.particle_id3, d.particle_id4)] + for d_id in dihedrals_to_delete: + del self._instances["dihedral"][d_id] + if "dihedral" in self._instances and not self._instances["dihedral"]: + del self._instances["dihedral"] + 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 @@ -415,6 +440,17 @@ def _get_templates_df(self, pmb_type): "central_particle": tpl.central_particle, "side_particle2": tpl.side_particle2, "parameters": parameters}) + elif pmb_type == "dihedral": + terms = [{key: value.to_quantity(self._units) for key, value in term.items()} + for term in tpl.terms] + rows.append({"pmb_type": tpl.pmb_type, + "name": tpl.name, + "dihedral_type": tpl.dihedral_type, + "particle_name1": tpl.particle_name1, + "particle_name2": tpl.particle_name2, + "particle_name3": tpl.particle_name3, + "particle_name4": tpl.particle_name4, + "terms": terms}) else: # Generic representation for other types rows.append(tpl.dict()) @@ -491,6 +527,9 @@ def _register_instance(self, instance): elif isinstance(instance, AngleInstance): pmb_type = "angle" iid = instance.angle_id + elif isinstance(instance, DihedralInstance): + pmb_type = "dihedral" + iid = instance.dihedral_id elif isinstance(instance, HydrogelInstance): pmb_type = "hydrogel" iid = instance.assembly_id @@ -574,6 +613,8 @@ def _update_instance(self, instance_id, pmb_type, attribute, value): allowed = ["particle_id1","particle_id2","added_to_engine"] elif pmb_type =="angle": allowed = ["particle_id1","particle_id2","particle_id3","added_to_engine"] + elif pmb_type =="dihedral": + allowed = ["particle_id1","particle_id2","particle_id3","particle_id4","added_to_engine"] else: allowed = [None] # No attributes allowed for other types if attribute not in allowed: @@ -747,6 +788,7 @@ def delete_instance(self, pmb_type, instance_id): if pmb_type == "particle": self._delete_bonds_of_particle(instance_id) self._delete_angles_of_particle(instance_id) + self._delete_dihedrals_of_particle(instance_id) # =============== FINAL DELETION STEP ====================== del self._instances[pmb_type][instance_id] if not self._instances[pmb_type]: @@ -826,6 +868,10 @@ def delete_template(self, pmb_type, name): if pmb_type == "angle": if name in self.espresso_angle_instances.keys(): del self.espresso_angle_instances[name] + # if it is a dihedral template delete also stored espresso dihedral instances + if pmb_type == "dihedral": + if name in self.espresso_dihedral_instances.keys(): + del self.espresso_dihedral_instances[name] # Delete empty groups if not self._templates[pmb_type]: del self._templates[pmb_type] diff --git a/pyMBE/storage/templates/dihedral.py b/pyMBE/storage/templates/dihedral.py new file mode 100644 index 0000000..1e2e228 --- /dev/null +++ b/pyMBE/storage/templates/dihedral.py @@ -0,0 +1,142 @@ +# +# 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, List, Literal, Optional +from ..base_type import PMBBaseModel +from ..pint_quantity import PintQuantity +from pydantic import Field + +class DihedralTemplate(PMBBaseModel): + """ + Template defining a dihedral (torsion) potential in the pyMBE database. + + Attributes: + pmb_type ('Literal["dihedral"]'): + Fixed type identifier for this template. Always "dihedral". + + name ('str'): + Unique name of the dihedral template, e.g., "A-B-C-D" where the + quadruplet is canonicalised under full reversal (see `make_dihedral_key`). + + dihedral_type ('str'): + Type of dihedral potential. Supported: "periodic" (the ESPResSo + ``Dihedral`` bond ``V(phi) = k*[1 - cos(n*phi - phi_0)]``). "tabulated" + is reserved for force fields that need Ryckaert-Bellemans or improper + dihedrals and is not implemented yet. + + particle_name1 ('Optional[str]'): + Name of the first particle in the dihedral quadruplet. + + particle_name2 ('Optional[str]'): + Name of the second particle in the dihedral quadruplet. + + particle_name3 ('Optional[str]'): + Name of the third particle in the dihedral quadruplet. + + particle_name4 ('Optional[str]'): + Name of the fourth particle in the dihedral quadruplet. + + terms ('List[Dict[str, PintQuantity]]'): + List of dihedral terms. GROMACS-style force fields stack several + ``[dihedrals]`` lines (different multiplicities) on the same + quadruplet, so a dihedral template holds a list of terms. For + ``dihedral_type == "periodic"`` each term is a dictionary with keys + ``k`` (energy), ``phi_0`` (dimensionless, radians) and ``mult`` + (dimensionless, the integer multiplicity). + + Notes: + - Values of the parameters are stored as PintQuantity objects for + unit-aware calculations. The multiplicity ``mult`` is stored as a + dimensionless PintQuantity so that all terms round-trip through the + same serialisation path. + - The canonical name is invariant under full reversal of the + quadruplet: ``make_dihedral_key("A", "B", "C", "D")`` equals + ``make_dihedral_key("D", "C", "B", "A")``. + - GROMACS proper dihedral type 1, ``V(phi) = k*[1 + cos(n*phi - phi_s)]``, + maps onto the ``"periodic"`` type with ``k`` unchanged, ``mult = n`` + and ``phi_0 = phi_s + pi`` (``phi_s`` converted from degrees to radians). + """ + pmb_type: Literal["dihedral"] = "dihedral" + name: str = Field(default="default") + dihedral_type: str + particle_name1: Optional[str] = None + particle_name2: Optional[str] = None + particle_name3: Optional[str] = None + particle_name4: Optional[str] = None + terms: List[Dict[str, PintQuantity]] + + @classmethod + def make_dihedral_key(cls, name1, name2, name3, name4): + """Return a canonical name for a dihedral between four particle names. + + A dihedral ``(a, b, c, d)`` describes the same torsion as its full + reversal ``(d, c, b, a)``. The canonical name is the lexicographically + smaller of the two orderings, joined by dashes. + + Args: + name1 ('str'): + Name of the first particle. + + name2 ('str'): + Name of the second particle. + + name3 ('str'): + Name of the third particle. + + name4 ('str'): + Name of the fourth particle. + + Returns: + ('str'): + Canonical dihedral name, e.g. "A-B-C-D". + """ + forward = (name1, name2, name3, name4) + reverse = (name4, name3, name2, name1) + n1, n2, n3, n4 = min(forward, reverse) + return f"{n1}-{n2}-{n3}-{n4}" + + def _make_name(self): + """Create canonical name using particle names.""" + if (not self.particle_name1 or not self.particle_name2 + or not self.particle_name3 or not self.particle_name4): + raise RuntimeError("Cannot generate dihedral name: one or more of " + "particle_name1..4 is missing.") + self.name = self.make_dihedral_key(self.particle_name1, + self.particle_name2, + self.particle_name3, + self.particle_name4) + + def get_terms(self, ureg): + """ + Retrieve the dihedral terms as dictionaries of Pint `Quantity` objects. + + Args: + ureg ('pint.UnitRegistry'): + Pint unit registry used to reconstruct physical quantities from storage. + + Returns: + 'List[Dict[str, pint.Quantity]]': + A list with one dictionary per term, mapping parameter names to + their corresponding unit-aware Pint quantities. + """ + pint_terms = [] + for term in self.terms: + pint_terms.append({key: value.to_quantity(ureg) + for key, value in term.items()}) + return pint_terms diff --git a/samples/setup_dihedral_potential.py b/samples/setup_dihedral_potential.py new file mode 100644 index 0000000..79a8057 --- /dev/null +++ b/samples/setup_dihedral_potential.py @@ -0,0 +1,180 @@ +# +# 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 dihedral (torsion) potential support in pyMBE. +# +# This script shows: +# 1. create_molecule with gen_dihedral=True +# - builds a linear chain of residues; dihedrals are auto-generated +# across the full molecule from the bond topology +# 2. a multi-term dihedral on a single quadruplet +# - GROMACS-style force fields stack several [dihedrals] lines (with +# different multiplicities) on the same four atoms; pyMBE registers one +# ESPResSo Dihedral object per term +# 3. validation of the ESPResSo <-> GROMACS proper-dihedral (type 1) mapping +# - ESPResSo Dihedral: V(phi) = k*[1 - cos(n*phi - phi_0)] +# - GROMACS type 1: V(phi) = k*[1 + cos(n*phi - phi_s)] +# - identical for phi_0 = phi_s + pi (phi_s converted from degrees) + + +import numpy as np +import espressomd +import pyMBE + +# ---------------------------------------------------------------------- +# Set up the ESPResSo system and pyMBE +# ---------------------------------------------------------------------- +box_l = [20] * 3 +espresso_system = espressomd.System(box_l=box_l) +pmb = pyMBE.pymbe_library(seed=42) +pmb.set_simulation_engine(espresso_system) + +# ---------------------------------------------------------------------- +# Particle and default bond definitions +# ---------------------------------------------------------------------- +for particle_name in ("A", "B", "C", "D"): + pmb.define_particle(name=particle_name, + z=0, + sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + +pmb.define_default_bond(bond_type='harmonic', + bond_parameters={'r_0': 0.4 * pmb.units.nm, + 'k': 40 * pmb.units('reduced_energy / reduced_length**2')}) + + +def dihedral_terms_on_particle(particle_id): + """Return the ESPResSo Dihedral interaction objects bonded to `particle_id`.""" + terms = [] + for bond_tuple in espresso_system.part.by_id(particle_id).bonds: + if len(bond_tuple) == 4 and type(bond_tuple[0]).__name__ == "Dihedral": + terms.append(bond_tuple[0]) + return terms + + +# ====================================================================== +# Demo 1: create_molecule with gen_dihedral=True +# ====================================================================== +# Topology: a linear chain A-B-C-D. The single bonded quadruplet +# (A, B, C, D) picks up the periodic dihedral defined below. +pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters={'k': 5 * pmb.units('reduced_energy'), + 'phi_0': np.pi * pmb.units(''), + 'mult': 1}, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + +pmb.define_residue(name="Res_A", central_bead="A", side_chains=[]) +pmb.define_residue(name="Res_B", central_bead="B", side_chains=[]) +pmb.define_residue(name="Res_C", central_bead="C", side_chains=[]) +pmb.define_residue(name="Res_D", central_bead="D", side_chains=[]) +pmb.define_molecule(name="Mol_1", residue_list=["Res_A", "Res_B", "Res_C", "Res_D"]) + +molecule_ids = pmb.create_molecule(name="Mol_1", + number_of_molecules=1, + box_l=box_l, + backbone_vector=[1.0, 0.0, 0.0], + use_default_bond=True, + gen_dihedral=True) +pmb.add_instances_to_engine() + +molecule_pids = pmb.db._find_instance_ids_by_attribute(pmb_type="particle", + attribute="molecule_id", + value=molecule_ids[0]) +print("\n############ Demo 1: linear molecule with gen_dihedral=True ############") +print(pmb.get_instances_df(pmb_type="dihedral")) +assert len(pmb.get_instances_df(pmb_type="dihedral")) == 1, "Demo 1: expected one dihedral" +assert len(dihedral_terms_on_particle(min(molecule_pids))) == 1 +print("Demo 1: one backbone dihedral generated OK") + +pmb.delete_instances_in_system(instance_id=molecule_ids[0], pmb_type="molecule") +pmb.db.delete_templates(pmb_type="dihedral") + + +# ====================================================================== +# Demo 2: multi-term dihedral on one quadruplet +# ====================================================================== +multi_term_parameters = [ + {'k': 4 * pmb.units('reduced_energy'), 'phi_0': 0.0 * pmb.units(''), 'mult': 1}, + {'k': 2 * pmb.units('reduced_energy'), 'phi_0': np.pi * pmb.units(''), 'mult': 3}, +] +pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=multi_term_parameters, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + +pids = [pmb.create_particle(name=name, box_l=box_l, number_of_particles=1)[0] + for name in ("A", "B", "C", "D")] +pmb.create_bond(particle_id1=pids[0], particle_id2=pids[1], use_default_bond=True) +pmb.create_bond(particle_id1=pids[1], particle_id2=pids[2], use_default_bond=True) +pmb.create_bond(particle_id1=pids[2], particle_id2=pids[3], use_default_bond=True) +pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3]) +pmb.add_instances_to_engine() + +print("\n############ Demo 2: multi-term dihedral ############") +print(f"registered {len(dihedral_terms_on_particle(pids[0]))} ESPResSo Dihedral objects " + f"for a {len(multi_term_parameters)}-term dihedral") +assert len(dihedral_terms_on_particle(pids[0])) == len(multi_term_parameters) +print("Demo 2: one Dihedral object per term OK") + +for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, pmb_type="particle") +pmb.db.delete_templates(pmb_type="dihedral") + + +# ====================================================================== +# Demo 3: ESPResSo <-> GROMACS proper dihedral (type 1) mapping +# ====================================================================== +# A GROMACS proper dihedral line "phi_s k n" (phi_s in degrees) maps onto a +# pyMBE "periodic" dihedral with the same k, mult = n and phi_0 = radians(phi_s) + pi. +gromacs_phi_s_deg = 180.0 +gromacs_k = 3.5 # reduced units here; kJ/mol in a real GROMACS topology +gromacs_n = 2 +phi_0 = np.deg2rad(gromacs_phi_s_deg) + np.pi + +pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters={'k': gromacs_k * pmb.units('reduced_energy'), + 'phi_0': phi_0 * pmb.units(''), + 'mult': gromacs_n}, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + +pids = [pmb.create_particle(name=name, box_l=box_l, number_of_particles=1)[0] + for name in ("A", "B", "C", "D")] +pmb.create_bond(particle_id1=pids[0], particle_id2=pids[1], use_default_bond=True) +pmb.create_bond(particle_id1=pids[1], particle_id2=pids[2], use_default_bond=True) +pmb.create_bond(particle_id1=pids[2], particle_id2=pids[3], use_default_bond=True) +pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3]) +pmb.add_instances_to_engine() + +params = dihedral_terms_on_particle(pids[0])[0].get_params() +print("\n############ Demo 3: GROMACS proper dihedral (type 1) mapping ############") +print(f" bend = {params['bend']} (GROMACS k = {gromacs_k})") +print(f" mult = {params['mult']} (GROMACS n = {gromacs_n})") +print(f" phase = {params['phase']:.6f} (radians(phi_s) + pi = {phi_0:.6f})") +assert np.isclose(params['bend'], gromacs_k) +assert params['mult'] == gromacs_n +assert np.isclose(params['phase'], phi_0) +print("Demo 3: GROMACS proper dihedral (type 1) mapping verified OK") + +for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, pmb_type="particle") +pmb.db.delete_templates(pmb_type="dihedral") diff --git a/testsuite/CTestTestfile.cmake b/testsuite/CTestTestfile.cmake index 687344a..ae657bb 100644 --- a/testsuite/CTestTestfile.cmake +++ b/testsuite/CTestTestfile.cmake @@ -59,6 +59,7 @@ 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 dihedral_tests.py) pymbe_add_test(PATH generate_perpendicular_vectors_test.py) pymbe_add_test(PATH define_and_create_molecules_unit_tests.py) pymbe_add_test(PATH create_molecule_position_test.py) diff --git a/testsuite/database_unit_tests.py b/testsuite/database_unit_tests.py index 23d3041..2965a22 100644 --- a/testsuite/database_unit_tests.py +++ b/testsuite/database_unit_tests.py @@ -27,9 +27,11 @@ 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.dihedral import DihedralInstance 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.dihedral import DihedralTemplate from pyMBE.storage.templates.hydrogel import HydrogelNode from pyMBE.storage.pint_quantity import PintQuantity from pyMBE.storage.reactions.reaction import Reaction, ReactionParticipant @@ -444,7 +446,16 @@ def test_instance_id_validators(self): self.assertRaises(ValueError, AngleInstance, **inputs) - + inputs = {"name":"A", + "dihedral_id":-1, + "particle_id1":1, + "particle_id2":2, + "particle_id3":3, + "particle_id4":4} + self.assertRaises(ValueError, + DihedralInstance, + **inputs) + def test_make_name_bond_template(self): inputs = {"bond_type": "harmonic", "parameters": {"r": PintQuantity(magnitude=1, @@ -465,7 +476,29 @@ def test_make_name_angle_template(self): angle_tpl = AngleTemplate(**inputs) self.assertRaises(RuntimeError, angle_tpl._make_name) - + + def test_make_name_dihedral_template(self): + inputs = {"dihedral_type": "periodic", + "terms": [{"k": PintQuantity(magnitude=1, + units="reduced_energy", + dimension="energy"), + "phi_0": PintQuantity(magnitude=1, + units="dimensionless", + dimension="dimensionless"), + "mult": PintQuantity(magnitude=1, + units="dimensionless", + dimension="dimensionless")}]} + dihedral_tpl = DihedralTemplate(**inputs) + self.assertRaises(RuntimeError, + dihedral_tpl._make_name) + + def test_canonical_dihedral_key(self): + """The canonical dihedral key is invariant under full reversal.""" + self.assertEqual(DihedralTemplate.make_dihedral_key("A", "B", "C", "D"), + DihedralTemplate.make_dihedral_key("D", "C", "B", "A")) + self.assertEqual(DihedralTemplate.make_dihedral_key("D", "C", "B", "A"), + "A-B-C-D") + def test_exceptions_pint_quantity(self): units = pint.UnitRegistry() inputs = {"q":1, diff --git a/testsuite/dihedral_tests.py b/testsuite/dihedral_tests.py new file mode 100644 index 0000000..17d58d2 --- /dev/null +++ b/testsuite/dihedral_tests.py @@ -0,0 +1,523 @@ +# +# 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 +box_l=[10]*3 +espresso_system = espressomd.System(box_l=box_l) + + +class Test(ut.TestCase): + def setUp(self): + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + + def define_templates(self, pmb): + for name in ('A', 'B', 'C', 'D'): + pmb.define_particle(name=name, + z=0, + sigma=0.4*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + + self.periodic_dihedral_params = { + 'k': 5 * pmb.units('reduced_energy'), + 'phi_0': np.pi * pmb.units(''), + 'mult': 1, + } + + self.periodic_dihedral_params_alt = { + 'k': 3 * pmb.units('reduced_energy'), + 'phi_0': np.pi / 2 * pmb.units(''), + 'mult': 2, + } + + self.multi_term_dihedral_params = [ + {'k': 4 * pmb.units('reduced_energy'), 'phi_0': 0.0 * pmb.units(''), 'mult': 1}, + {'k': 2 * pmb.units('reduced_energy'), 'phi_0': np.pi * pmb.units(''), 'mult': 3}, + ] + + self.harmonic_bond_params = { + 'r_0': 0.4 * pmb.units.nm, + 'k': 400 * pmb.units('reduced_energy / reduced_length**2'), + } + + def get_dihedral_objects(self, particle_id1): + """ + Returns the list of dihedral interaction objects stored in ESPResSo for a + given first particle. Dihedral bonds are added to the first particle in + ESPResSo and carry three partner particles (tuple arity 4). + """ + bonds = espresso_system.part.by_id(particle_id1).bonds + dihedral_objects = [] + for bond_tuple in bonds: + if len(bond_tuple) == 4: + dihedral_objects.append(bond_tuple[0]) + return dihedral_objects + + def check_dihedral_setup(self, dihedral_objects, input_terms): + """ + Checks that pyMBE sets up the dihedral interaction objects correctly. + + Args: + dihedral_objects: list of ESPResSo dihedral interaction objects. + input_terms('list[dict]'): list of the term dictionaries passed to pyMBE. + """ + self.assertEqual(len(dihedral_objects), len(input_terms), + msg="pyMBE did not register one ESPResSo Dihedral object per term") + for dihedral_object in dihedral_objects: + self.assertIn("dihedral", str(type(dihedral_object)).lower(), + msg="pyMBE does not store the correct type of dihedral interaction") + actual = sorted((d.get_params()['bend'], + d.get_params()['mult'], + round(d.get_params()['phase'], 10)) for d in dihedral_objects) + expected = sorted((term['k'].m_as('reduced_energy'), + int(term['mult']), + round(float(term['phi_0'].magnitude), 10)) for term in input_terms) + for (b_a, m_a, p_a), (b_e, m_e, p_e) in zip(actual, expected): + np.testing.assert_almost_equal(actual=b_a, desired=b_e) + self.assertEqual(m_a, m_e) + np.testing.assert_almost_equal(actual=p_a, desired=p_e) + + def build_linear_chain(self, pmb): + """Create particles A-B-C-D and the three backbone bonds.""" + pids = [] + for name in ('A', 'B', 'C', 'D'): + pids.append(pmb.create_particle(name=name, + box_l=box_l, + number_of_particles=1)[0]) + pmb.create_bond(particle_id1=pids[0], particle_id2=pids[1]) + pmb.create_bond(particle_id1=pids[1], particle_id2=pids[2]) + pmb.create_bond(particle_id1=pids[2], particle_id2=pids[3]) + return pids + + def test_dihedral_setup_periodic(self): + """ + Unit test to check the setup of a periodic dihedral potential 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'], ['C', 'D']]) + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + + pids = self.build_linear_chain(pmb) + pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3]) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + + dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + self.assertEqual(len(dihedral_objects), 1, "No dihedral object found on the first particle") + self.check_dihedral_setup(dihedral_objects=dihedral_objects, + input_terms=[self.periodic_dihedral_params]) + + for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + def test_multi_term_dihedral(self): + """ + Unit test to check that a multi-term dihedral registers one ESPResSo + Dihedral object per term, all attached to the same quadruplet. + """ + 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'], ['C', 'D']]) + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.multi_term_dihedral_params, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + + pids = self.build_linear_chain(pmb) + pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3]) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + + dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + self.assertEqual(len(dihedral_objects), 2, "Expected two Dihedral objects for a two-term dihedral") + self.check_dihedral_setup(dihedral_objects=dihedral_objects, + input_terms=self.multi_term_dihedral_params) + self.assertEqual(len(pmb.get_instances_df(pmb_type="dihedral")), 1) + + for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + def test_dihedral_instance_reuses_cached_object(self): + """ + Test that cached ESPResSo dihedral interactions are reused for the same template. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + dihedral_template = pmb.get_dihedral_template(name1="A", name2="B", name3="C", name4="D") + + pmb.set_simulation_engine(espresso_system) + first = pmb.simulation_engine._get_dihedral_instance(dihedral_template=dihedral_template) + second = pmb.simulation_engine._get_dihedral_instance(dihedral_template=dihedral_template) + + self.assertIs(first, second) + self.assertEqual(len(pmb.db.espresso_dihedral_instances), 1) + + pmb.db.delete_templates(pmb_type="dihedral") + + def test_default_dihedral(self): + """ + Unit test to check the setup of a default dihedral potential. + """ + 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'], ['C', 'D']]) + pmb.define_default_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params) + + pids = self.build_linear_chain(pmb) + pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3], + use_default_dihedral=True) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + + dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + self.assertEqual(len(dihedral_objects), 1) + self.check_dihedral_setup(dihedral_objects=dihedral_objects, + input_terms=[self.periodic_dihedral_params]) + + for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + def test_specific_dihedral_over_default(self): + """ + Test that a specific dihedral template is preferred over the default. + """ + 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'], ['C', 'D']]) + pmb.define_default_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params_alt) + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + + pids = self.build_linear_chain(pmb) + pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3], + use_default_dihedral=True) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + + dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + self.check_dihedral_setup(dihedral_objects=dihedral_objects, + input_terms=[self.periodic_dihedral_params]) + + for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + def test_dihedral_raised_exceptions(self): + """ + Unit test to check that dihedral-related methods raise the correct exceptions. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + for callback in [pmb.define_dihedral_potential, pmb.define_default_dihedral_potential]: + with self.subTest(msg=f'using method {callback.__qualname__}()'): + self.check_dihedral_exceptions(callback, pmb) + + def check_dihedral_exceptions(self, callback, pmb): + good_term = {'k': 5 * pmb.units('reduced_energy'), + 'phi_0': np.pi * pmb.units(''), + 'mult': 1} + + def call(dihedral_type, dihedral_parameters): + kwargs = {"dihedral_type": dihedral_type, + "dihedral_parameters": dihedral_parameters} + if callback == pmb.define_dihedral_potential: + kwargs["particle_quadruplets"] = [('A', 'B', 'C', 'D')] + return kwargs + + # unknown dihedral type + np.testing.assert_raises(NotImplementedError, callback, + **call('quartic', good_term)) + # tabulated is reserved but not implemented yet + np.testing.assert_raises(NotImplementedError, callback, + **call('tabulated', good_term)) + # missing k + np.testing.assert_raises(ValueError, callback, + **call('periodic', {'phi_0': np.pi * pmb.units(''), 'mult': 1})) + # missing phi_0 + np.testing.assert_raises(ValueError, callback, + **call('periodic', {'k': 5 * pmb.units('reduced_energy'), 'mult': 1})) + # missing mult + np.testing.assert_raises(ValueError, callback, + **call('periodic', {'k': 5 * pmb.units('reduced_energy'), + 'phi_0': np.pi * pmb.units('')})) + # duplicate quadruplet in define_dihedral_potential + if callback == pmb.define_dihedral_potential: + np.testing.assert_raises(RuntimeError, + pmb.define_dihedral_potential, + dihedral_type="periodic", + dihedral_parameters=good_term, + particle_quadruplets=[('A', 'B', 'C', 'D'), + ('D', 'C', 'B', 'A')]) + + def test_sanity_get_dihedral_template(self): + """ + Tests that get_dihedral_template raises ValueError when no template is found. + """ + pmb = pyMBE.pymbe_library(seed=51) + np.testing.assert_raises(ValueError, + pmb.get_dihedral_template, + name1="W", name2="X", name3="Y", name4="Z") + + def test_canonical_dihedral_name(self): + """ + Test that dihedral names are canonical (invariant under full reversal). + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params, + particle_quadruplets=[('D', 'C', 'B', 'A')]) + + tpl = pmb.get_dihedral_template(name1="A", name2="B", name3="C", name4="D") + self.assertEqual(tpl.name, "A-B-C-D") + tpl2 = pmb.get_dihedral_template(name1="D", name2="C", name3="B", name4="A") + self.assertEqual(tpl2.name, "A-B-C-D") + + pmb.db.delete_templates(pmb_type="dihedral") + + def test_dihedral_without_bonds_raises_error(self): + """ + Test that creating a dihedral without bonds between the particles raises a ValueError. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + + pids = [] + for name in ('A', 'B', 'C', 'D'): + pids.append(pmb.create_particle(name=name, + box_l=box_l, + number_of_particles=1)[0]) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + with self.assertRaises(ValueError): + pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3]) + + for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle") + pmb.db.delete_templates(pmb_type="dihedral") + + def test_dihedral_with_partial_bonds_raises_error(self): + """ + Test that creating a dihedral with a missing bond in the i-j-k-l chain raises a ValueError. + """ + 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_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + + pids = [] + for name in ('A', 'B', 'C', 'D'): + pids.append(pmb.create_particle(name=name, + box_l=box_l, + number_of_particles=1)[0]) + pmb.create_bond(particle_id1=pids[0], particle_id2=pids[1]) + pmb.create_bond(particle_id1=pids[1], particle_id2=pids[2]) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + with self.assertRaises(ValueError): + pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3]) + + for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + def test_generate_dihedrals_for_entity_without_particles_returns(self): + """ + Test that auto-generating dihedrals for an empty entity returns without changes. + """ + pmb = pyMBE.pymbe_library(seed=42) + self.define_templates(pmb) + + pmb._generate_dihedrals_for_entity(entity_id=999, + entity_id_col="residue_id") + self.assertEqual(len(pmb.get_instances_df(pmb_type="dihedral")), 0) + + def test_generate_dihedrals_for_entity_creates_dihedral_from_bonds(self): + """ + Test that auto-generating dihedrals creates one dihedral from a bonded quadruplet. + """ + 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'], ['C', 'D']]) + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + + pids = self.build_linear_chain(pmb) + for particle_id in pids: + pmb.db._update_instance(instance_id=particle_id, + pmb_type="particle", + attribute="residue_id", + value=0) + + pmb._generate_dihedrals_for_entity(entity_id=0, + entity_id_col="residue_id") + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + + dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + self.assertEqual(len(dihedral_objects), 1) + self.assertEqual(len(pmb.get_instances_df(pmb_type="dihedral")), 1) + + for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + def test_generate_dihedrals_for_entity_skips_missing_templates(self): + """ + Test that auto-generating dihedrals skips bonded quadruplets without a matching 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'], ['C', 'D']]) + + pids = self.build_linear_chain(pmb) + for particle_id in pids: + pmb.db._update_instance(instance_id=particle_id, + pmb_type="particle", + attribute="residue_id", + value=0) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + pmb._generate_dihedrals_for_entity(entity_id=0, + entity_id_col="residue_id") + + self.assertEqual(len(self.get_dihedral_objects(particle_id1=pids[0])), 0) + self.assertEqual(len(pmb.get_instances_df(pmb_type="dihedral")), 0) + + for inst_id in pids: + pmb.delete_instances_in_system(instance_id=inst_id, + pmb_type="particle") + pmb.db.delete_templates(pmb_type="bond") + + def test_create_molecule_with_gen_dihedral_generates_dihedrals(self): + """ + Test that create_molecule generates backbone dihedrals automatically when gen_dihedral 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'], ['C', 'D']]) + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params, + particle_quadruplets=[('A', 'B', 'C', 'D')]) + 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_residue(name="R_D", central_bead="D", side_chains=[]) + pmb.define_molecule(name="M_dihedral", + residue_list=["R_A", "R_B", "R_C", "R_D"]) + + molecule_ids = pmb.create_molecule(name="M_dihedral", + number_of_molecules=1, + box_l=box_l, + backbone_vector=[1.0, 0.0, 0.0], + gen_dihedral=True) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + + self.assertEqual(len(pmb.get_instances_df(pmb_type="dihedral")), 1) + + pmb.delete_instances_in_system(instance_id=molecule_ids[0], + pmb_type="molecule") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + +if __name__ == '__main__': + ut.main() diff --git a/testsuite/test_io_database.py b/testsuite/test_io_database.py index 7667f36..66528ab 100644 --- a/testsuite/test_io_database.py +++ b/testsuite/test_io_database.py @@ -320,6 +320,72 @@ def test_io_angle_templates(self): pd.testing.assert_frame_equal(pmb.get_templates_df(pmb_type="angle"), new_pmb.get_templates_df(pmb_type="angle")) + def test_io_dihedral_templates(self): + """ + Checks that information in the pyMBE database about + dihedral templates is stored to file and loaded back properly, + including multi-term dihedrals. + """ + pmb = pyMBE.pymbe_library(seed=42) + parameters1 = {"k": 5.0 * pmb.units.reduced_energy, + "phi_0": 1.0 * pmb.units(""), + "mult": 1} + parameters2 = [{"k": 3.0 * pmb.units.reduced_energy, "phi_0": 0.0 * pmb.units(""), "mult": 1}, + {"k": 2.0 * pmb.units.reduced_energy, "phi_0": 0.5 * pmb.units(""), "mult": 3}] + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=parameters1, + particle_quadruplets=[("Z", "X", "Z", "Y"), + ("X", "Z", "Y", "W")]) + pmb.define_default_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=parameters2) + new_pmb = pyMBE.pymbe_library(23) + with tempfile.TemporaryDirectory() as tmp_directory: + pmb.save_database(tmp_directory) + new_pmb.load_database(tmp_directory) + pd.testing.assert_frame_equal(pmb.get_templates_df(pmb_type="dihedral"), + new_pmb.get_templates_df(pmb_type="dihedral")) + + def test_io_dihedral_instances(self): + """ + Checks that dihedral instances are stored to file and loaded back properly. + """ + pmb = pyMBE.pymbe_library(seed=42) + for name in ("A", "B", "C", "D"): + pmb.define_particle(name=name, + sigma=3.5 * pmb.units.reduced_length, + cutoff=4 * pmb.units.reduced_length, + offset=0 * pmb.units.reduced_length, + epsilon=0.2 * pmb.units.reduced_energy, + z=0) + parameters = {"k": 100.0 * pmb.units.reduced_energy / (pmb.units.reduced_length**2), + "r_0": 1.0 * pmb.units.reduced_length} + pmb.define_bond(bond_type="harmonic", + bond_parameters=parameters, + particle_pairs=[["A", "B"], ["B", "C"], ["C", "D"]]) + pmb.define_dihedral_potential(dihedral_type="periodic", + dihedral_parameters={"k": 5.0 * pmb.units.reduced_energy, + "phi_0": 1.0 * pmb.units(""), + "mult": 2}, + particle_quadruplets=[("A", "B", "C", "D")]) + pids = [pmb.create_particle(name=name, box_l=box_l, number_of_particles=1)[0] + for name in ("A", "B", "C", "D")] + pmb.create_bond(particle_id1=pids[0], particle_id2=pids[1]) + pmb.create_bond(particle_id1=pids[1], particle_id2=pids[2]) + pmb.create_bond(particle_id1=pids[2], particle_id2=pids[3]) + pmb.create_dihedral_potential(particle_id1=pids[0], + particle_id2=pids[1], + particle_id3=pids[2], + particle_id4=pids[3]) + + new_pmb = pyMBE.pymbe_library(23) + with tempfile.TemporaryDirectory() as tmp_directory: + pmb.save_database(tmp_directory) + new_pmb.load_database(tmp_directory) + pd.testing.assert_frame_equal(pmb.get_instances_df(pmb_type="dihedral"), + new_pmb.get_instances_df(pmb_type="dihedral")) + pd.testing.assert_frame_equal(pmb.get_templates_df(pmb_type="dihedral"), + new_pmb.get_templates_df(pmb_type="dihedral")) + def test_io_residues_templates(self): """ Checks that information in the pyMBE database about From c3b10a2636ea040d47c14604a8f58a7702f422e9 Mon Sep 17 00:00:00 2001 From: Raj Kumar Rajaram Baskaran <12559850+rajkumarrb@users.noreply.github.com> Date: Wed, 2 Sep 2026 20:52:07 +0100 Subject: [PATCH 2/4] Fix simulation_builder id-collision test for dihedral lookup add_instances_to_engine now performs a fifth _find_instance_ids_by_attribute call for pmb_type='dihedral', mirroring the bond and angle blocks. The pre-existing test_espresso_engine_add_instances_id_collision mocked that method with a four-element side_effect list, so the dihedral lookup raised StopIteration. Add the missing entry. --- testsuite/simulation_builder_tests.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/testsuite/simulation_builder_tests.py b/testsuite/simulation_builder_tests.py index 5a4d577..73cfa38 100644 --- a/testsuite/simulation_builder_tests.py +++ b/testsuite/simulation_builder_tests.py @@ -154,7 +154,8 @@ def test_espresso_engine_add_instances_id_collision(self): local_pmb.set_simulation_engine(espresso_system) engine = local_pmb.simulation_engine engine.db = Mock() - engine.db._find_instance_ids_by_attribute.side_effect = [[1], [2], [], []] + # particle (missing), particle (added), bond, angle, dihedral lookups + engine.db._find_instance_ids_by_attribute.side_effect = [[1], [2], [], [], []] engine._get_particle_ids_in_espresso = Mock(return_value=[1, 2]) engine.update_particle_id = Mock() engine._add_particle = Mock() From f3c8b86f2a8f5ec5d442f12e9e83766714599db2 Mon Sep 17 00:00:00 2001 From: Raj Kumar Rajaram Baskaran <12559850+rajkumarrb@users.noreply.github.com> Date: Wed, 28 Oct 2026 08:33:40 +0000 Subject: [PATCH 3/4] Use GROMACS dihedral convention in the API and fix ESPResSo bond attachment Dihedral parameters (k, phi_s, mult) now follow the IUPAC/GROMACS type 1 convention, V(phi) = k*[1 + cos(n*phi - phi_s)] with phi = 0 at the cis conformation. The ESPResSo engine converts them internally to the ESPResSo Dihedral form (phase = -(phi_s + pi)) and documents the conversion. Attach the ESPResSo Dihedral bond to the second particle of the quadruplet, as required by ESPResSo, instead of the first one, which made ESPResSo evaluate a different torsion. Add tests that compare energies and forces set up in ESPResSo with the analytical expression: energy and forces vs. finite differences, momentum conservation, NVE energy conservation, relaxation to the minimum, GROMACS type 1 parameters, reversed quadruplets, periodic boundaries, and generation of dihedrals for residues, molecules and hydrogels. --- CHANGELOG.md | 3 +- pyMBE/pyMBE.py | 26 +- pyMBE/simulation_builder/espresso_engine.py | 22 +- pyMBE/storage/templates/dihedral.py | 13 +- samples/setup_dihedral_potential.py | 36 +- testsuite/database_unit_tests.py | 2 +- testsuite/dihedral_tests.py | 448 +++++++++++++++++++- testsuite/test_io_database.py | 8 +- 8 files changed, 495 insertions(+), 63 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 019daec..2979280 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,7 +14,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Function `create_dihedral_potential` to create instances of dihedral potentials in pyMBE. (#156) - Function `get_dihedral_template` to retrieve a template of a dihedral potential in the pyMBE database. (#156) - `gen_dihedral` keyword argument in `create_residue`, `create_molecule` and `create_hydrogel` to auto-generate dihedral potentials from the bond topology. (#156) -- Sample script `samples/setup_dihedral_potential.py` and test suite `testsuite/dihedral_tests.py` for the new functionality. (#156) +- Dihedral potentials follow the IUPAC/GROMACS convention at the API level: users pass the parameters of the GROMACS proper dihedral (type 1), `V(phi) = k*[1 + cos(n*phi - phi_s)]` with `phi = 0` at the cis conformation, as `k`, `phi_s` and `mult`. The conversion to the ESPResSo `Dihedral` interaction is applied internally by the ESPResSo engine adapter. (#156) +- Sample script `samples/setup_dihedral_potential.py` and test suite `testsuite/dihedral_tests.py` for the new functionality, including tests that compare the energies and forces set up in ESPResSo with the analytical GROMACS expression. (#156) - Simulation-engine abstraction layer to decouple pyMBE molecular-building logic from simulation-engine-specific implementation details. (#151) - Added `EspressoEngine` to delegate ESPResSo-specific operations, including reaction methods, setup utilities, and particle/bond/angle insertion. (#151) - Added `add_instances_to_engine()` to transfer pyMBE particles, bonds, and angles to the selected simulation engine. (#151) diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py index e1d582c..d8d7f85 100644 --- a/pyMBE/pyMBE.py +++ b/pyMBE/pyMBE.py @@ -1662,7 +1662,7 @@ def _build_dihedral_terms(self, dihedral_type, dihedral_parameters): dictionary is treated as a one-term dihedral. For the "periodic" type each term must contain: - "k" ('pint.Quantity'): Torsion stiffness with dimensions of energy. - - "phi_0" ('pint.Quantity'): Equilibrium phase, dimensionless (radians). + - "phi_s" ('pint.Quantity'): Phase of the GROMACS proper dihedral (type 1), dimensionless (radians). The dihedral angle follows the IUPAC convention (cis = 0, trans = pi). - "mult" ('int'): Multiplicity (number of minima). Returns: @@ -1685,8 +1685,8 @@ def _build_dihedral_terms(self, dihedral_type, dihedral_parameters): for term in term_list: if "k" not in term: raise ValueError("Magnitude of the dihedral potential (k) is missing") - if "phi_0" not in term: - raise ValueError("Equilibrium phase of the dihedral potential (phi_0) is missing") + if "phi_s" not in term: + raise ValueError("Phase of the dihedral potential (phi_s) is missing") if "mult" not in term: raise ValueError("Multiplicity of the dihedral potential (mult) is missing") mult = term["mult"] @@ -1695,7 +1695,7 @@ def _build_dihedral_terms(self, dihedral_type, dihedral_parameters): terms.append({"k": PintQuantity.from_quantity(q=term["k"], expected_dimension="energy", ureg=self.units), - "phi_0": PintQuantity.from_quantity(q=term["phi_0"], + "phi_s": PintQuantity.from_quantity(q=term["phi_s"], expected_dimension="dimensionless", ureg=self.units), "mult": PintQuantity.from_quantity(q=mult, @@ -1709,8 +1709,8 @@ def define_dihedral_potential(self, dihedral_type, dihedral_parameters, particle Args: dihedral_type ('str'): - Type of dihedral potential. Supported: "periodic", which maps onto - the ESPResSo ``Dihedral`` bond ``V(phi) = k*[1 - cos(n*phi - phi_0)]``. + Type of dihedral potential. Supported: "periodic", the GROMACS proper + dihedral (type 1) ``V(phi) = k*[1 + cos(n*phi - phi_s)]``. dihedral_parameters ('dict' or 'list[dict]'): Parameters of the dihedral potential, or a list of such @@ -1718,16 +1718,20 @@ def define_dihedral_potential(self, dihedral_type, dihedral_parameters, particle on the same quadruplet. For the "periodic" type each term must contain: - "k" ('pint.Quantity'): Torsion stiffness with dimensions of energy. - - "phi_0" ('pint.Quantity'): Equilibrium phase, dimensionless (radians). + - "phi_s" ('pint.Quantity'): Phase of the GROMACS proper dihedral (type 1), dimensionless (radians). The dihedral angle follows the IUPAC convention (cis = 0, trans = pi). - "mult" ('int'): Multiplicity (number of minima). particle_quadruplets ('list[tuple[str,str,str,str]]'): List of (particle1, particle2, particle3, particle4) name quadruplets. Notes: - GROMACS proper dihedral type 1, ``V(phi) = k*[1 + cos(n*phi - phi_s)]``, - is recovered by passing ``k`` unchanged, ``mult = n`` and - ``phi_0 = phi_s + pi`` (with ``phi_s`` converted from degrees to radians). + Parameters follow the GROMACS proper dihedral (type 1) convention, + ``V(phi) = k*[1 + cos(n*phi - phi_s)]``: pass ``k``, ``mult = n`` and + ``phi_s`` (converted from degrees to radians) exactly as written in a + GROMACS topology. The dihedral angle ``phi`` follows the IUPAC/GROMACS + convention (cis = 0, trans = pi, right-handed rotation positive). + Any conversion needed by a simulation engine is applied internally by + the engine adapter, see `EspressoSimulation._create_dihedral_instance`. """ terms = self._build_dihedral_terms(dihedral_type=dihedral_type, dihedral_parameters=dihedral_parameters) @@ -1754,7 +1758,7 @@ def define_default_dihedral_potential(self, dihedral_type, dihedral_parameters): Type of dihedral potential. Supported: "periodic". dihedral_parameters ('dict' or 'list[dict]'): - Parameters of the dihedral potential (k, phi_0, mult), or a list + Parameters of the dihedral potential (k, phi_s, mult), or a list of such dictionaries for a multi-term dihedral. """ terms = self._build_dihedral_terms(dihedral_type=dihedral_type, diff --git a/pyMBE/simulation_builder/espresso_engine.py b/pyMBE/simulation_builder/espresso_engine.py index dbae054..d60243d 100644 --- a/pyMBE/simulation_builder/espresso_engine.py +++ b/pyMBE/simulation_builder/espresso_engine.py @@ -73,8 +73,9 @@ def _add_dihedral(self,particle_id1,particle_id2,particle_id3,particle_id4, dihe pmb_type="dihedral") espresso_dihedral_insts=self._get_dihedral_instance(dihedral_template=dihedral_tpl) + # ESPResSo stores a dihedral bond on the second particle of the quadruplet: p2.add_bond((dihedral, p1, p3, p4)) for espresso_dihedral_inst in espresso_dihedral_insts: - self.espresso_system.part.by_id(particle_id1).add_bond((espresso_dihedral_inst, particle_id2, particle_id3, particle_id4)) + self.espresso_system.part.by_id(particle_id2).add_bond((espresso_dihedral_inst, particle_id1, particle_id3, particle_id4)) self.db._update_instance(instance_id=dihedral_inst.dihedral_id, pmb_type='dihedral', attribute='added_to_engine', @@ -301,10 +302,18 @@ def _create_dihedral_instance(self, dihedral_type, dihedral_terms): Args: dihedral_type ('str'): Type of dihedral potential ("periodic"). - dihedral_terms ('list[dict]'): List of term dictionaries (k, phi_0, mult). + dihedral_terms ('list[dict]'): List of term dictionaries (k, phi_s, mult). Returns: ('list[espressomd.interactions.BondedInteraction]'): The ESPResSo dihedral interaction objects. + + Notes: + pyMBE stores dihedrals in the IUPAC/GROMACS convention, ``V = k*[1 + cos(n*phi - phi_s)]`` + with phi = 0 at the cis conformation. ESPResSo implements ``V = k*[1 - cos(n*phi_E - phase)]`` + and, for a bond stored on the second particle of the quadruplet, its angle ``phi_E`` has the + opposite sign to the IUPAC angle. The conversion applied here is therefore + ``phase = -(phi_s + pi)``: the ``pi`` turns ``1 + cos`` into ``1 - cos`` and the minus + sign compensates the sign of ``phi_E``. """ from espressomd import interactions @@ -312,9 +321,14 @@ def _create_dihedral_instance(self, dihedral_type, dihedral_terms): dihedral_insts = [] for term in dihedral_terms: k = term["k"].m_as("reduced_energy") - phi_0 = float(term["phi_0"].magnitude) + phi_s = float(term["phi_s"].magnitude) mult = int(round(float(term["mult"].magnitude))) - dihedral_insts.append(interactions.Dihedral(bend=k, mult=mult, phase=phi_0)) + # pyMBE uses the GROMACS proper dihedral (type 1), V = k[1 + cos(n*phi - phi_s)], + # with the IUPAC angle phi (cis = 0). ESPResSo implements V = k[1 - cos(n*phi_E - phase)]: + # - 1 + cos(x) = 1 - cos(x - pi), so the phase is shifted by pi: phi_0 = phi_s + pi + # - with the bond stored on the second particle (see _add_dihedral), ESPResSo's angle + # phi_E has the opposite sign to the IUPAC angle, so phase = -phi_0 + dihedral_insts.append(interactions.Dihedral(bend=k, mult=mult, phase=-(phi_s + np.pi))) return dihedral_insts if dihedral_type == "tabulated": raise NotImplementedError("tabulated dihedrals not yet supported") diff --git a/pyMBE/storage/templates/dihedral.py b/pyMBE/storage/templates/dihedral.py index 1e2e228..bf7a0ef 100644 --- a/pyMBE/storage/templates/dihedral.py +++ b/pyMBE/storage/templates/dihedral.py @@ -35,8 +35,8 @@ class DihedralTemplate(PMBBaseModel): quadruplet is canonicalised under full reversal (see `make_dihedral_key`). dihedral_type ('str'): - Type of dihedral potential. Supported: "periodic" (the ESPResSo - ``Dihedral`` bond ``V(phi) = k*[1 - cos(n*phi - phi_0)]``). "tabulated" + Type of dihedral potential. Supported: "periodic" (the GROMACS + proper dihedral, type 1, ``V(phi) = k*[1 + cos(n*phi - phi_s)]``). "tabulated" is reserved for force fields that need Ryckaert-Bellemans or improper dihedrals and is not implemented yet. @@ -57,7 +57,7 @@ class DihedralTemplate(PMBBaseModel): ``[dihedrals]`` lines (different multiplicities) on the same quadruplet, so a dihedral template holds a list of terms. For ``dihedral_type == "periodic"`` each term is a dictionary with keys - ``k`` (energy), ``phi_0`` (dimensionless, radians) and ``mult`` + ``k`` (energy), ``phi_s`` (dimensionless, radians) and ``mult`` (dimensionless, the integer multiplicity). Notes: @@ -68,9 +68,10 @@ class DihedralTemplate(PMBBaseModel): - The canonical name is invariant under full reversal of the quadruplet: ``make_dihedral_key("A", "B", "C", "D")`` equals ``make_dihedral_key("D", "C", "B", "A")``. - - GROMACS proper dihedral type 1, ``V(phi) = k*[1 + cos(n*phi - phi_s)]``, - maps onto the ``"periodic"`` type with ``k`` unchanged, ``mult = n`` - and ``phi_0 = phi_s + pi`` (``phi_s`` converted from degrees to radians). + - The parameters are stored in the GROMACS proper dihedral (type 1) + convention, ``V(phi) = k*[1 + cos(n*phi - phi_s)]`` with ``phi_s`` + in radians and ``phi`` in the IUPAC convention (cis = 0). The + simulation-engine adapter converts them to the engine's own form. """ pmb_type: Literal["dihedral"] = "dihedral" name: str = Field(default="default") diff --git a/samples/setup_dihedral_potential.py b/samples/setup_dihedral_potential.py index 79a8057..4ae389d 100644 --- a/samples/setup_dihedral_potential.py +++ b/samples/setup_dihedral_potential.py @@ -27,9 +27,10 @@ # different multiplicities) on the same four atoms; pyMBE registers one # ESPResSo Dihedral object per term # 3. validation of the ESPResSo <-> GROMACS proper-dihedral (type 1) mapping -# - ESPResSo Dihedral: V(phi) = k*[1 - cos(n*phi - phi_0)] -# - GROMACS type 1: V(phi) = k*[1 + cos(n*phi - phi_s)] -# - identical for phi_0 = phi_s + pi (phi_s converted from degrees) +# - pyMBE takes the GROMACS type 1 parameters as they appear in a topology: +# V(phi) = k*[1 + cos(n*phi - phi_s)], phi in the IUPAC convention (cis = 0) +# - the ESPResSo engine adapter converts them to the ESPResSo Dihedral form +# V(phi) = k*[1 - cos(n*phi - phase)] internally, phase = -(phi_s + pi) import numpy as np @@ -59,7 +60,7 @@ def dihedral_terms_on_particle(particle_id): - """Return the ESPResSo Dihedral interaction objects bonded to `particle_id`.""" + """Return the ESPResSo Dihedral interaction objects bonded to `particle_id` (the second particle of the quadruplet).""" terms = [] for bond_tuple in espresso_system.part.by_id(particle_id).bonds: if len(bond_tuple) == 4 and type(bond_tuple[0]).__name__ == "Dihedral": @@ -74,7 +75,7 @@ def dihedral_terms_on_particle(particle_id): # (A, B, C, D) picks up the periodic dihedral defined below. pmb.define_dihedral_potential(dihedral_type="periodic", dihedral_parameters={'k': 5 * pmb.units('reduced_energy'), - 'phi_0': np.pi * pmb.units(''), + 'phi_s': np.pi * pmb.units(''), 'mult': 1}, particle_quadruplets=[('A', 'B', 'C', 'D')]) @@ -98,7 +99,7 @@ def dihedral_terms_on_particle(particle_id): print("\n############ Demo 1: linear molecule with gen_dihedral=True ############") print(pmb.get_instances_df(pmb_type="dihedral")) assert len(pmb.get_instances_df(pmb_type="dihedral")) == 1, "Demo 1: expected one dihedral" -assert len(dihedral_terms_on_particle(min(molecule_pids))) == 1 +assert len(dihedral_terms_on_particle(sorted(molecule_pids)[1])) == 1 print("Demo 1: one backbone dihedral generated OK") pmb.delete_instances_in_system(instance_id=molecule_ids[0], pmb_type="molecule") @@ -109,8 +110,8 @@ def dihedral_terms_on_particle(particle_id): # Demo 2: multi-term dihedral on one quadruplet # ====================================================================== multi_term_parameters = [ - {'k': 4 * pmb.units('reduced_energy'), 'phi_0': 0.0 * pmb.units(''), 'mult': 1}, - {'k': 2 * pmb.units('reduced_energy'), 'phi_0': np.pi * pmb.units(''), 'mult': 3}, + {'k': 4 * pmb.units('reduced_energy'), 'phi_s': 0.0 * pmb.units(''), 'mult': 1}, + {'k': 2 * pmb.units('reduced_energy'), 'phi_s': np.pi * pmb.units(''), 'mult': 3}, ] pmb.define_dihedral_potential(dihedral_type="periodic", dihedral_parameters=multi_term_parameters, @@ -128,9 +129,9 @@ def dihedral_terms_on_particle(particle_id): pmb.add_instances_to_engine() print("\n############ Demo 2: multi-term dihedral ############") -print(f"registered {len(dihedral_terms_on_particle(pids[0]))} ESPResSo Dihedral objects " +print(f"registered {len(dihedral_terms_on_particle(pids[1]))} ESPResSo Dihedral objects " f"for a {len(multi_term_parameters)}-term dihedral") -assert len(dihedral_terms_on_particle(pids[0])) == len(multi_term_parameters) +assert len(dihedral_terms_on_particle(pids[1])) == len(multi_term_parameters) print("Demo 2: one Dihedral object per term OK") for inst_id in pids: @@ -141,16 +142,17 @@ def dihedral_terms_on_particle(particle_id): # ====================================================================== # Demo 3: ESPResSo <-> GROMACS proper dihedral (type 1) mapping # ====================================================================== -# A GROMACS proper dihedral line "phi_s k n" (phi_s in degrees) maps onto a -# pyMBE "periodic" dihedral with the same k, mult = n and phi_0 = radians(phi_s) + pi. +# A GROMACS proper dihedral line "phi_s k n" (phi_s in degrees) is passed to pyMBE as it is: +# same k, mult = n and phi_s converted to radians. The ESPResSo engine applies the conversion +# to the ESPResSo Dihedral form internally (phase = -(phi_s + pi)). gromacs_phi_s_deg = 180.0 gromacs_k = 3.5 # reduced units here; kJ/mol in a real GROMACS topology gromacs_n = 2 -phi_0 = np.deg2rad(gromacs_phi_s_deg) + np.pi +phi_s = np.deg2rad(gromacs_phi_s_deg) pmb.define_dihedral_potential(dihedral_type="periodic", dihedral_parameters={'k': gromacs_k * pmb.units('reduced_energy'), - 'phi_0': phi_0 * pmb.units(''), + 'phi_s': phi_s * pmb.units(''), 'mult': gromacs_n}, particle_quadruplets=[('A', 'B', 'C', 'D')]) @@ -165,14 +167,14 @@ def dihedral_terms_on_particle(particle_id): particle_id4=pids[3]) pmb.add_instances_to_engine() -params = dihedral_terms_on_particle(pids[0])[0].get_params() +params = dihedral_terms_on_particle(pids[1])[0].get_params() print("\n############ Demo 3: GROMACS proper dihedral (type 1) mapping ############") print(f" bend = {params['bend']} (GROMACS k = {gromacs_k})") print(f" mult = {params['mult']} (GROMACS n = {gromacs_n})") -print(f" phase = {params['phase']:.6f} (radians(phi_s) + pi = {phi_0:.6f})") +print(f" phase = {params['phase']:.6f} (-(radians(phi_s) + pi) = {-(phi_s + np.pi):.6f}, converted by the ESPResSo engine)") assert np.isclose(params['bend'], gromacs_k) assert params['mult'] == gromacs_n -assert np.isclose(params['phase'], phi_0) +assert np.isclose(params["phase"], -(phi_s + np.pi)) print("Demo 3: GROMACS proper dihedral (type 1) mapping verified OK") for inst_id in pids: diff --git a/testsuite/database_unit_tests.py b/testsuite/database_unit_tests.py index 2965a22..68854f0 100644 --- a/testsuite/database_unit_tests.py +++ b/testsuite/database_unit_tests.py @@ -482,7 +482,7 @@ def test_make_name_dihedral_template(self): "terms": [{"k": PintQuantity(magnitude=1, units="reduced_energy", dimension="energy"), - "phi_0": PintQuantity(magnitude=1, + "phi_s": PintQuantity(magnitude=1, units="dimensionless", dimension="dimensionless"), "mult": PintQuantity(magnitude=1, diff --git a/testsuite/dihedral_tests.py b/testsuite/dihedral_tests.py index 17d58d2..832aa14 100644 --- a/testsuite/dihedral_tests.py +++ b/testsuite/dihedral_tests.py @@ -18,13 +18,17 @@ # Import pyMBE and other libraries import pyMBE +from pyMBE.lib.lattice import DiamondLattice import numpy as np import unittest as ut import espressomd +import itertools # Create an instance of the ESPResSo system box_l=[10]*3 espresso_system = espressomd.System(box_l=box_l) +espresso_system.time_step = 0.002 +espresso_system.cell_system.skin = 0.4 class Test(ut.TestCase): @@ -41,19 +45,19 @@ def define_templates(self, pmb): self.periodic_dihedral_params = { 'k': 5 * pmb.units('reduced_energy'), - 'phi_0': np.pi * pmb.units(''), + 'phi_s': np.pi * pmb.units(''), 'mult': 1, } self.periodic_dihedral_params_alt = { 'k': 3 * pmb.units('reduced_energy'), - 'phi_0': np.pi / 2 * pmb.units(''), + 'phi_s': np.pi / 2 * pmb.units(''), 'mult': 2, } self.multi_term_dihedral_params = [ - {'k': 4 * pmb.units('reduced_energy'), 'phi_0': 0.0 * pmb.units(''), 'mult': 1}, - {'k': 2 * pmb.units('reduced_energy'), 'phi_0': np.pi * pmb.units(''), 'mult': 3}, + {'k': 4 * pmb.units('reduced_energy'), 'phi_s': 0.0 * pmb.units(''), 'mult': 1}, + {'k': 2 * pmb.units('reduced_energy'), 'phi_s': np.pi * pmb.units(''), 'mult': 3}, ] self.harmonic_bond_params = { @@ -61,13 +65,13 @@ def define_templates(self, pmb): 'k': 400 * pmb.units('reduced_energy / reduced_length**2'), } - def get_dihedral_objects(self, particle_id1): + def get_dihedral_objects(self, particle_id2): """ Returns the list of dihedral interaction objects stored in ESPResSo for a - given first particle. Dihedral bonds are added to the first particle in - ESPResSo and carry three partner particles (tuple arity 4). + given second particle. Dihedral bonds are added to the second particle of + the quadruplet in ESPResSo and carry three partner particles (tuple arity 4). """ - bonds = espresso_system.part.by_id(particle_id1).bonds + bonds = espresso_system.part.by_id(particle_id2).bonds dihedral_objects = [] for bond_tuple in bonds: if len(bond_tuple) == 4: @@ -77,6 +81,7 @@ def get_dihedral_objects(self, particle_id1): def check_dihedral_setup(self, dihedral_objects, input_terms): """ Checks that pyMBE sets up the dihedral interaction objects correctly. + pyMBE takes GROMACS type 1 parameters; the ESPResSo phase is -(phi_s + pi). Args: dihedral_objects: list of ESPResSo dihedral interaction objects. @@ -92,7 +97,7 @@ def check_dihedral_setup(self, dihedral_objects, input_terms): round(d.get_params()['phase'], 10)) for d in dihedral_objects) expected = sorted((term['k'].m_as('reduced_energy'), int(term['mult']), - round(float(term['phi_0'].magnitude), 10)) for term in input_terms) + round(-(float(term["phi_s"].magnitude) + np.pi), 10)) for term in input_terms) for (b_a, m_a, p_a), (b_e, m_e, p_e) in zip(actual, expected): np.testing.assert_almost_equal(actual=b_a, desired=b_e) self.assertEqual(m_a, m_e) @@ -132,7 +137,7 @@ def test_dihedral_setup_periodic(self): pmb.set_simulation_engine(espresso_system) pmb.add_instances_to_engine() - dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + dihedral_objects = self.get_dihedral_objects(particle_id2=pids[1]) self.assertEqual(len(dihedral_objects), 1, "No dihedral object found on the first particle") self.check_dihedral_setup(dihedral_objects=dihedral_objects, input_terms=[self.periodic_dihedral_params]) @@ -166,7 +171,7 @@ def test_multi_term_dihedral(self): pmb.set_simulation_engine(espresso_system) pmb.add_instances_to_engine() - dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + dihedral_objects = self.get_dihedral_objects(particle_id2=pids[1]) self.assertEqual(len(dihedral_objects), 2, "Expected two Dihedral objects for a two-term dihedral") self.check_dihedral_setup(dihedral_objects=dihedral_objects, input_terms=self.multi_term_dihedral_params) @@ -221,7 +226,7 @@ def test_default_dihedral(self): pmb.set_simulation_engine(espresso_system) pmb.add_instances_to_engine() - dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + dihedral_objects = self.get_dihedral_objects(particle_id2=pids[1]) self.assertEqual(len(dihedral_objects), 1) self.check_dihedral_setup(dihedral_objects=dihedral_objects, input_terms=[self.periodic_dihedral_params]) @@ -257,7 +262,7 @@ def test_specific_dihedral_over_default(self): pmb.set_simulation_engine(espresso_system) pmb.add_instances_to_engine() - dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + dihedral_objects = self.get_dihedral_objects(particle_id2=pids[1]) self.check_dihedral_setup(dihedral_objects=dihedral_objects, input_terms=[self.periodic_dihedral_params]) @@ -280,7 +285,7 @@ def test_dihedral_raised_exceptions(self): def check_dihedral_exceptions(self, callback, pmb): good_term = {'k': 5 * pmb.units('reduced_energy'), - 'phi_0': np.pi * pmb.units(''), + 'phi_s': np.pi * pmb.units(''), 'mult': 1} def call(dihedral_type, dihedral_parameters): @@ -298,14 +303,14 @@ def call(dihedral_type, dihedral_parameters): **call('tabulated', good_term)) # missing k np.testing.assert_raises(ValueError, callback, - **call('periodic', {'phi_0': np.pi * pmb.units(''), 'mult': 1})) - # missing phi_0 + **call('periodic', {'phi_s': np.pi * pmb.units(''), 'mult': 1})) + # missing phi_s np.testing.assert_raises(ValueError, callback, **call('periodic', {'k': 5 * pmb.units('reduced_energy'), 'mult': 1})) # missing mult np.testing.assert_raises(ValueError, callback, **call('periodic', {'k': 5 * pmb.units('reduced_energy'), - 'phi_0': np.pi * pmb.units('')})) + 'phi_s': np.pi * pmb.units('')})) # duplicate quadruplet in define_dihedral_potential if callback == pmb.define_dihedral_potential: np.testing.assert_raises(RuntimeError, @@ -443,7 +448,7 @@ def test_generate_dihedrals_for_entity_creates_dihedral_from_bonds(self): pmb.set_simulation_engine(espresso_system) pmb.add_instances_to_engine() - dihedral_objects = self.get_dihedral_objects(particle_id1=pids[0]) + dihedral_objects = self.get_dihedral_objects(particle_id2=pids[1]) self.assertEqual(len(dihedral_objects), 1) self.assertEqual(len(pmb.get_instances_df(pmb_type="dihedral")), 1) @@ -475,7 +480,7 @@ def test_generate_dihedrals_for_entity_skips_missing_templates(self): pmb._generate_dihedrals_for_entity(entity_id=0, entity_id_col="residue_id") - self.assertEqual(len(self.get_dihedral_objects(particle_id1=pids[0])), 0) + self.assertEqual(len(self.get_dihedral_objects(particle_id2=pids[1])), 0) self.assertEqual(len(pmb.get_instances_df(pmb_type="dihedral")), 0) for inst_id in pids: @@ -518,6 +523,411 @@ def test_create_molecule_with_gen_dihedral_generates_dihedrals(self): pmb.db.delete_templates(pmb_type="dihedral") pmb.db.delete_templates(pmb_type="bond") + def test_create_nested_residue_with_gen_dihedral_generates_dihedrals(self): + """ + Test that create_residue generates the dihedral of a chain A-B-C-D formed by nested residues. + """ + 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'], ['C', 'D']]) + pmb.define_default_dihedral_potential(dihedral_type="periodic", + dihedral_parameters=self.periodic_dihedral_params) + pmb.define_residue(name="R3", central_bead="C", side_chains=["D"]) + pmb.define_residue(name="R2", central_bead="B", side_chains=["R3"]) + pmb.define_residue(name="R1", central_bead="A", side_chains=["R2"]) + + residue_ids = pmb.create_residue(name="R1", + box_l=box_l, + gen_dihedral=True) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + + self.assertEqual(len(pmb.get_instances_df(pmb_type="dihedral")), 1) + n_engine_dihedrals = sum(1 for particle in espresso_system.part.all() + for bond in particle.bonds if len(bond) == 4) + self.assertEqual(n_engine_dihedrals, 1) + + pmb.delete_instances_in_system(instance_id=residue_ids, + pmb_type="residue") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + def test_create_hydrogel_with_gen_dihedral_generates_chain_dihedrals(self): + """ + Test that create_hydrogel generates the dihedrals internal to every chain, i.e. mpc-3 per chain, + and none involving the crosslinker nodes. + """ + mpc = 5 + pmb = pyMBE.pymbe_library(seed=42) + bond_length = 0.355 * pmb.units.nm + pmb.define_particle(name="N", sigma=bond_length, epsilon=1 * pmb.units("reduced_energy")) + pmb.define_particle(name="C", sigma=bond_length, epsilon=1 * pmb.units("reduced_energy")) + pmb.define_residue(name="Res", central_bead="C", side_chains=[]) + pmb.define_molecule(name="chain", residue_list=["Res"] * mpc) + pmb.define_bond(bond_type="harmonic", + bond_parameters={"r_0": bond_length, + "k": 400 * pmb.units("reduced_energy / reduced_length**2")}, + particle_pairs=[["C", "C"], ["C", "N"]]) + pmb.define_default_dihedral_potential(dihedral_type="periodic", + dihedral_parameters={'k': 5 * pmb.units('reduced_energy'), + 'phi_s': np.pi * pmb.units(''), + 'mult': 1}) + lattice = DiamondLattice(mpc, bond_length) + hydrogel_box_l = [lattice.box_l] * 3 + try: + espresso_system.box_l = hydrogel_box_l + pmb.initialize_lattice_builder(lattice) + node_topology = [{"particle_name": "N", "lattice_index": index} for index in lattice.indices] + chain_topology = [{"node_start": str(lattice.indices[connectivity[0]]), + "node_end": str(lattice.indices[connectivity[1]]), + "molecule_name": "chain"} for connectivity in lattice.connectivity] + pmb.define_hydrogel("gel", node_topology, chain_topology) + pmb.create_hydrogel("gel", hydrogel_box_l, gen_dihedral=True) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + + dihedrals = pmb.db.get_instances(pmb_type="dihedral") + particle_names = {pid: inst.name for pid, inst in pmb.db.get_instances(pmb_type="particle").items()} + self.assertEqual(len(dihedrals), len(lattice.connectivity) * (mpc - 3)) + for dihedral in dihedrals.values(): + names = [particle_names[getattr(dihedral, f"particle_id{i}")] for i in (1, 2, 3, 4)] + self.assertNotIn("N", names, msg="crosslinker-adjacent dihedrals are not generated") + n_engine_dihedrals = sum(1 for particle in espresso_system.part.all() + for bond in particle.bonds if len(bond) == 4) + self.assertEqual(n_engine_dihedrals, len(dihedrals)) + finally: + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + espresso_system.box_l = box_l + + +def iupac_dihedral(p1, p2, p3, p4): + """ + Independent reference for the dihedral angle of the quadruplet p1-p2-p3-p4 + in the IUPAC/GROMACS convention (cis = 0, trans = pi, right-handed positive). + """ + b1, b2, b3 = p2 - p1, p3 - p2, p4 - p3 + n1, n2 = np.cross(b1, b2), np.cross(b2, b3) + m1 = np.cross(n1, b2 / np.linalg.norm(b2)) + return np.arctan2(np.dot(m1, n2), np.dot(n1, n2)) + + +def periodic_energy(positions, terms): + """Analytic GROMACS type 1 energy sum_t k_t [1 + cos(n_t phi - phi_s_t)] of a quadruplet, terms = [(k, n, phi_s), ...].""" + phi = iupac_dihedral(*positions) + return sum(k * (1 + np.cos(n * phi - phi_s)) for k, n, phi_s in terms) + + +class DihedralPhysicsTest(ut.TestCase): + """ + Checks that the dihedral potential set up by pyMBE in ESPResSo is the right + physical potential, using references computed independently with numpy: + energies, forces, momentum conservation, energy conservation in MD, + relaxation to the minimum, symmetry, GROMACS mapping and periodic images. + """ + K = 5.0 + TERMS = [(1, 0.0), (1, 1.0), (2, 0.3), (3, 2.0), (3, np.pi)] # (mult, phi_s) + + def setUp(self): + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + espresso_system.integrator.set_vv() + espresso_system.thermostat.turn_off() + self.rng = np.random.default_rng(42) + + def build(self, terms): + """ + Builds A-B-C-D with harmonic backbone bonds and a dihedral with the given + terms [(k, n, phi_s), ...] and returns (pmb, pids, r_0 in reduced units). + """ + pmb = pyMBE.pymbe_library(seed=42) + for name in ('A', 'B', 'C', 'D'): + pmb.define_particle(name=name, z=0, sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + r_0 = 0.4 * pmb.units.nm + pmb.define_bond(bond_type="harmonic", + bond_parameters={'r_0': r_0, + 'k': 400 * pmb.units('reduced_energy / reduced_length**2')}, + particle_pairs=[['A', 'B'], ['B', 'C'], ['C', 'D']]) + pmb.define_dihedral_potential( + dihedral_type="periodic", + dihedral_parameters=[{'k': k * pmb.units('reduced_energy'), + 'phi_s': phi_s * pmb.units(''), + 'mult': n} for k, n, phi_s in terms], + particle_quadruplets=[('A', 'B', 'C', 'D')]) + pids = [pmb.create_particle(name=name, box_l=box_l, number_of_particles=1)[0] + for name in ('A', 'B', 'C', 'D')] + for i in range(3): + pmb.create_bond(particle_id1=pids[i], particle_id2=pids[i + 1]) + pmb.create_dihedral_potential(particle_id1=pids[0], particle_id2=pids[1], + particle_id3=pids[2], particle_id4=pids[3]) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + return pmb, pids, r_0.m_as('reduced_length') + + @staticmethod + def geometry(r_0, theta, alpha=1.7, beta=1.9, center=(5., 5., 5.)): + """ + Positions of A, B, C, D with all three bond lengths exactly r_0 (so the harmonic + bonds contribute no energy or force). alpha/beta are the polar angles of the + A and D bonds relative to the B->C axis, theta is the rotation about it. + """ + B = np.array(center, dtype=float) + C = B + r_0 * np.array([0., 0., 1.]) + A = B + r_0 * np.array([np.sin(alpha), 0., np.cos(alpha)]) + D = C + r_0 * np.array([np.sin(beta) * np.cos(theta), np.sin(beta) * np.sin(theta), np.cos(beta)]) + return np.array([A, B, C, D]) + + @staticmethod + def place(pids, positions): + for pid, pos in zip(pids, positions): + espresso_system.part.by_id(pid).pos = pos + + @staticmethod + def cleanup(pmb, pids): + for pid in pids: + pmb.delete_instances_in_system(instance_id=pid, pmb_type="particle") + pmb.db.delete_templates(pmb_type="dihedral") + pmb.db.delete_templates(pmb_type="bond") + + def energy_and_forces(self, pids): + espresso_system.integrator.run(0, recalc_forces=True) + energy = espresso_system.analysis.energy()['bonded'] + forces = np.array([espresso_system.part.by_id(pid).f for pid in pids]) + return energy, forces + + def test_energy_matches_analytic_formula(self): + """ESPResSo energy equals the GROMACS expression k[1 + cos(n*phi - phi_s)] for several (n, phi_s) and geometries.""" + for n, phi_s in self.TERMS: + pmb, pids, r_0 = self.build([(self.K, n, phi_s)]) + for theta in np.linspace(-2.9, 2.9, 11): + pos = self.geometry(r_0, theta) + self.place(pids, pos) + energy, _ = self.energy_and_forces(pids) + np.testing.assert_allclose(energy, periodic_energy(pos, [(self.K, n, phi_s)]), + atol=1e-8, err_msg=f"n={n}, phi_s={phi_s}, theta={theta}") + self.cleanup(pmb, pids) + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + + def test_energy_is_not_symmetric_in_sign_of_phi(self): + """For phi_s not in {0, pi}, phi and -phi differ in energy: guards the sign convention.""" + phi_s = 1.0 + pmb, pids, r_0 = self.build([(self.K, 1, phi_s)]) + energies = [] + for theta in (0.8, -0.8): + pos = self.geometry(r_0, theta) + self.place(pids, pos) + energies.append(self.energy_and_forces(pids)[0]) + np.testing.assert_allclose(energies[-1], periodic_energy(pos, [(self.K, 1, phi_s)]), atol=1e-8) + self.assertGreater(abs(energies[0] - energies[1]), 0.1) + self.cleanup(pmb, pids) + + def test_multi_term_energy_is_the_sum_of_terms(self): + """A stacked multi-term dihedral adds the energies of its terms.""" + terms = [(4.0, 1, 0.0), (2.0, 3, np.pi), (1.5, 2, 0.7)] + pmb, pids, r_0 = self.build(terms) + for theta in np.linspace(-2.8, 2.8, 9): + pos = self.geometry(r_0, theta) + self.place(pids, pos) + np.testing.assert_allclose(self.energy_and_forces(pids)[0], periodic_energy(pos, terms), atol=1e-8) + self.cleanup(pmb, pids) + + def test_forces_match_numerical_gradient(self): + """ESPResSo forces equal -grad of the analytic energy (central finite differences).""" + terms = [(self.K, 2, 0.3), (1.5, 1, 1.0)] + pmb, pids, r_0 = self.build(terms) + h = 1e-6 + for theta in (-2.3, -1.1, 0.7, 1.9, 2.6): + pos = self.geometry(r_0, theta, alpha=1.5, beta=2.0) + self.place(pids, pos) + _, forces = self.energy_and_forces(pids) + reference = np.zeros_like(pos) + for i in range(4): + for d in range(3): + plus, minus = pos.copy(), pos.copy() + plus[i, d] += h + minus[i, d] -= h + reference[i, d] = -(periodic_energy(plus, terms) - periodic_energy(minus, terms)) / (2 * h) + np.testing.assert_allclose(forces, reference, atol=1e-5, err_msg=f"theta={theta}") + self.cleanup(pmb, pids) + + def test_forces_conserve_momentum_and_angular_momentum(self): + """Net force and net torque of an internal torsion vanish.""" + pmb, pids, r_0 = self.build([(self.K, 3, 0.4)]) + for theta in (-2.0, 0.6, 1.4): + pos = self.geometry(r_0, theta, alpha=1.3, beta=2.2) + self.place(pids, pos) + _, forces = self.energy_and_forces(pids) + self.assertGreater(np.abs(forces).max(), 1e-3, "test needs non-zero forces to be meaningful") + np.testing.assert_allclose(forces.sum(axis=0), 0., atol=1e-9) + np.testing.assert_allclose(np.cross(pos, forces).sum(axis=0), 0., atol=1e-9) + self.cleanup(pmb, pids) + + def test_nve_energy_conservation(self): + """Velocity-Verlet NVE dynamics conserves total energy (forces consistent with the energy).""" + pmb, pids, r_0 = self.build([(self.K, 2, 0.3), (2.0, 1, 1.0)]) + self.place(pids, self.geometry(r_0, 1.2, alpha=1.7, beta=1.9)) + for pid in pids: + espresso_system.part.by_id(pid).v = 0.3 * self.rng.normal(size=3) + espresso_system.integrator.run(0, recalc_forces=True) + e_0 = espresso_system.analysis.energy()['total'] + drift = 0. + for _ in range(40): + espresso_system.integrator.run(100) + drift = max(drift, abs(espresso_system.analysis.energy()['total'] - e_0)) + self.assertLess(drift, 5e-3 * max(1., abs(e_0)), msg=f"energy drift {drift} from E_0={e_0}") + self.cleanup(pmb, pids) + + def test_relaxation_reaches_the_minimum_with_the_right_sign(self): + """ + Steepest descent from a twisted conformation ends at the minimum of the GROMACS + type 1 potential, phi_min = phi_s + pi for n = 1 (IUPAC angle), not at -phi_min. + This depends on the sign convention and on the forces, not on the energy formula. + """ + for phi_min in (2.0, -2.0): + phi_s = phi_min - np.pi if phi_min > 0 else phi_min + np.pi + pmb, pids, r_0 = self.build([(self.K, 1, phi_s)]) + start = 1.0 if phi_min > 0 else -1.0 + self.place(pids, self.geometry(r_0, start)) + espresso_system.integrator.set_steepest_descent(f_max=1e-4, gamma=1e-3, max_displacement=0.01) # gamma*k_bond << 2 for stability + espresso_system.integrator.run(5000) + espresso_system.integrator.set_vv() + pos = np.array([espresso_system.part.by_id(pid).pos for pid in pids]) + # unfold the chain around B (bond lengths are far below half the box) + pos = pos[1] + (pos - pos[1] + np.array(box_l) / 2) % np.array(box_l) - np.array(box_l) / 2 + phi = iupac_dihedral(*pos) + self.assertAlmostEqual(phi, phi_min, places=3) + self.assertLess(espresso_system.analysis.energy()['bonded'], 1e-6) + self.cleanup(pmb, pids) + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + + def test_reversed_quadruplet_gives_the_same_physics(self): + """Creating the dihedral on (D, C, B, A) hits the same template and gives the same energy and forces.""" + terms = [(self.K, 2, 0.3), (1.5, 1, 1.0)] + pmb, pids, r_0 = self.build(terms) + pos = self.geometry(r_0, 0.9, alpha=1.4, beta=2.1) + self.place(pids, pos) + e_forward, f_forward = self.energy_and_forces(pids) + self.cleanup(pmb, pids) + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + + pmb = pyMBE.pymbe_library(seed=42) + for name in ('A', 'B', 'C', 'D'): + pmb.define_particle(name=name, z=0, sigma=0.4 * pmb.units.nm, + epsilon=1 * pmb.units('reduced_energy')) + pmb.define_bond(bond_type="harmonic", + bond_parameters={'r_0': 0.4 * pmb.units.nm, + 'k': 400 * pmb.units('reduced_energy / reduced_length**2')}, + particle_pairs=[['A', 'B'], ['B', 'C'], ['C', 'D']]) + pmb.define_dihedral_potential( + dihedral_type="periodic", + dihedral_parameters=[{'k': k * pmb.units('reduced_energy'), 'phi_s': phi_s * pmb.units(''), 'mult': n} + for k, n, phi_s in terms], + particle_quadruplets=[('A', 'B', 'C', 'D')]) + pids = [pmb.create_particle(name=name, box_l=box_l, number_of_particles=1)[0] + for name in ('A', 'B', 'C', 'D')] + for i in range(3): + pmb.create_bond(particle_id1=pids[i], particle_id2=pids[i + 1]) + pmb.create_dihedral_potential(particle_id1=pids[3], particle_id2=pids[2], + particle_id3=pids[1], particle_id4=pids[0]) + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + self.place(pids, pos) + e_reverse, f_reverse = self.energy_and_forces(pids) + np.testing.assert_allclose(e_reverse, e_forward, atol=1e-10) + np.testing.assert_allclose(f_reverse, f_forward, atol=1e-9) + np.testing.assert_allclose(e_reverse, periodic_energy(pos, terms), atol=1e-8) + self.cleanup(pmb, pids) + + def test_gromacs_proper_dihedral_mapping(self): + """ + A GROMACS proper dihedral (type 1), V = k[1 + cos(n*phi - phi_s)], given to pyMBE exactly + as in a topology (k, mult = n, phi_s in radians) is reproduced by the ESPResSo engine. + """ + for n in (1, 2, 3): + for phi_s_deg in (0., 180., 60., -90.): + phi_s = np.radians(phi_s_deg) + pmb, pids, r_0 = self.build([(self.K, n, phi_s)]) + for theta in (-2.4, -0.9, 0.5, 1.6, 2.7): + pos = self.geometry(r_0, theta) + self.place(pids, pos) + gromacs = self.K * (1 + np.cos(n * iupac_dihedral(*pos) - phi_s)) + np.testing.assert_allclose(self.energy_and_forces(pids)[0], gromacs, atol=1e-8, + err_msg=f"n={n}, phi_s={phi_s_deg}, theta={theta}") + self.cleanup(pmb, pids) + espresso_system.part.clear() + espresso_system.bonded_inter.clear() + + def test_dihedral_across_periodic_boundary(self): + """A torsion whose particles are on opposite sides of the box has the same energy and forces.""" + terms = [(self.K, 2, 0.3)] + pmb, pids, r_0 = self.build(terms) + results = [] + for center in ((5., 5., 5.), (box_l[0] - 0.4 * r_0, 0.3 * r_0, box_l[2] - 0.2 * r_0)): + pos = self.geometry(r_0, 0.9, alpha=1.4, beta=2.1, center=center) + self.place(pids, pos) + if center[0] != 5.: + self.assertTrue(np.any(pos > box_l[0]) or np.any(pos < 0.), + "test geometry should straddle the periodic boundary") + results.append(self.energy_and_forces(pids)) + np.testing.assert_allclose(results[-1][0], periodic_energy(pos, terms), atol=1e-8) + np.testing.assert_allclose(results[1][0], results[0][0], atol=1e-8) + np.testing.assert_allclose(results[1][1], results[0][1], atol=1e-8) + self.cleanup(pmb, pids) + + def count_generated_dihedrals(self, n_particles, bonds): + """Builds a bonded graph with a default dihedral and returns the generated dihedral quadruplets.""" + pmb = pyMBE.pymbe_library(seed=42) + pmb.define_particle(name='A', z=0, sigma=0.4 * pmb.units.nm, epsilon=1 * pmb.units('reduced_energy')) + pmb.define_default_bond(bond_type="harmonic", + bond_parameters={'r_0': 0.4 * pmb.units.nm, + 'k': 400 * pmb.units('reduced_energy / reduced_length**2')}) + pmb.define_default_dihedral_potential(dihedral_type="periodic", + dihedral_parameters={'k': 5 * pmb.units('reduced_energy'), + 'phi_s': np.pi * pmb.units(''), + 'mult': 1}) + pids = pmb.create_particle(name='A', box_l=box_l, number_of_particles=n_particles) + for i, j in bonds: + pmb.create_bond(particle_id1=pids[i], particle_id2=pids[j], use_default_bond=True) + pmb.db._update_instance(instance_id=pids[i], pmb_type="particle", attribute="residue_id", value=0) + pmb.db._update_instance(instance_id=pids[j], pmb_type="particle", attribute="residue_id", value=0) + pmb._generate_dihedrals_for_entity(entity_id=0, entity_id_col="residue_id") + index = {pid: i for i, pid in enumerate(pids)} + pmb.set_simulation_engine(espresso_system) + pmb.add_instances_to_engine() + found = [tuple(index[getattr(d, f"particle_id{k}")] for k in (1, 2, 3, 4)) + for d in pmb.db.get_instances(pmb_type="dihedral").values()] + self.cleanup(pmb, pids) + return found + + @staticmethod + def reference_torsions(n_particles, bonds): + """Independent count: simple 4-vertex paths over the bond graph, a path and its reverse being one torsion.""" + edges = {frozenset(b) for b in bonds} + paths = {p for p in itertools.permutations(range(n_particles), 4) + if all(frozenset(p[i:i + 2]) in edges for i in range(3))} + return {min(p, p[::-1]) for p in paths} + + def test_gen_dihedral_matches_graph_enumeration(self): + """Auto-generation finds exactly the torsions of the bond graph: chain, branch, rings, K4.""" + topologies = {"linear chain of 6": (6, [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5)]), + "branched": (5, [(0, 1), (1, 2), (2, 3), (2, 4)]), + "4-ring": (4, [(0, 1), (1, 2), (2, 3), (3, 0)]), + "triangle (no torsion)": (3, [(0, 1), (1, 2), (2, 0)]), + "K4": (4, [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)])} + for label, (n_particles, bonds) in topologies.items(): + found = self.count_generated_dihedrals(n_particles, bonds) + reference = self.reference_torsions(n_particles, bonds) + self.assertEqual(len(found), len(reference), msg=f"{label}: number of dihedrals") + self.assertEqual({min(q, q[::-1]) for q in found}, reference, msg=f"{label}: quadruplets") + espresso_system.part.clear() + if __name__ == '__main__': ut.main() diff --git a/testsuite/test_io_database.py b/testsuite/test_io_database.py index 66528ab..1768187 100644 --- a/testsuite/test_io_database.py +++ b/testsuite/test_io_database.py @@ -328,10 +328,10 @@ def test_io_dihedral_templates(self): """ pmb = pyMBE.pymbe_library(seed=42) parameters1 = {"k": 5.0 * pmb.units.reduced_energy, - "phi_0": 1.0 * pmb.units(""), + "phi_s": 1.0 * pmb.units(""), "mult": 1} - parameters2 = [{"k": 3.0 * pmb.units.reduced_energy, "phi_0": 0.0 * pmb.units(""), "mult": 1}, - {"k": 2.0 * pmb.units.reduced_energy, "phi_0": 0.5 * pmb.units(""), "mult": 3}] + parameters2 = [{"k": 3.0 * pmb.units.reduced_energy, "phi_s": 0.0 * pmb.units(""), "mult": 1}, + {"k": 2.0 * pmb.units.reduced_energy, "phi_s": 0.5 * pmb.units(""), "mult": 3}] pmb.define_dihedral_potential(dihedral_type="periodic", dihedral_parameters=parameters1, particle_quadruplets=[("Z", "X", "Z", "Y"), @@ -364,7 +364,7 @@ def test_io_dihedral_instances(self): particle_pairs=[["A", "B"], ["B", "C"], ["C", "D"]]) pmb.define_dihedral_potential(dihedral_type="periodic", dihedral_parameters={"k": 5.0 * pmb.units.reduced_energy, - "phi_0": 1.0 * pmb.units(""), + "phi_s": 1.0 * pmb.units(""), "mult": 2}, particle_quadruplets=[("A", "B", "C", "D")]) pids = [pmb.create_particle(name=name, box_l=box_l, number_of_particles=1)[0] From 791daeee9eb7da9f356f9d25bdac392e595ba93f Mon Sep 17 00:00:00 2001 From: Raj Kumar Rajaram Baskaran <12559850+rajkumarrb@users.noreply.github.com> Date: Wed, 28 Oct 2026 08:52:40 +0000 Subject: [PATCH 4/4] Reference pull request number in the changelog --- CHANGELOG.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2979280..1ecbb02 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,14 +8,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ## Added -- Support to setup dihedral (torsion) potentials with pyMBE, completing the bonded-potential support tracked in #143. All flexible pyMBE templates support dihedral potentials: molecules, peptides, residues (including residues with nested residues) and the internal chains of hydrogels. (#156) -- Template and instance `Dihedral` to store information about dihedral potentials in the pyMBE database, including support for multi-term (stacked multiplicity) dihedrals on the same quadruplet. (#156) -- Functions `define_dihedral_potential` and `define_default_dihedral_potential` to define templates of dihedral potentials in pyMBE. (#156) -- Function `create_dihedral_potential` to create instances of dihedral potentials in pyMBE. (#156) -- Function `get_dihedral_template` to retrieve a template of a dihedral potential in the pyMBE database. (#156) -- `gen_dihedral` keyword argument in `create_residue`, `create_molecule` and `create_hydrogel` to auto-generate dihedral potentials from the bond topology. (#156) -- Dihedral potentials follow the IUPAC/GROMACS convention at the API level: users pass the parameters of the GROMACS proper dihedral (type 1), `V(phi) = k*[1 + cos(n*phi - phi_s)]` with `phi = 0` at the cis conformation, as `k`, `phi_s` and `mult`. The conversion to the ESPResSo `Dihedral` interaction is applied internally by the ESPResSo engine adapter. (#156) -- Sample script `samples/setup_dihedral_potential.py` and test suite `testsuite/dihedral_tests.py` for the new functionality, including tests that compare the energies and forces set up in ESPResSo with the analytical GROMACS expression. (#156) +- Support to setup dihedral (torsion) potentials with pyMBE, completing the bonded-potential support tracked in #143. All flexible pyMBE templates support dihedral potentials: molecules, peptides, residues (including residues with nested residues) and the internal chains of hydrogels. (#157) +- Template and instance `Dihedral` to store information about dihedral potentials in the pyMBE database, including support for multi-term (stacked multiplicity) dihedrals on the same quadruplet. (#157) +- Functions `define_dihedral_potential` and `define_default_dihedral_potential` to define templates of dihedral potentials in pyMBE. (#157) +- Function `create_dihedral_potential` to create instances of dihedral potentials in pyMBE. (#157) +- Function `get_dihedral_template` to retrieve a template of a dihedral potential in the pyMBE database. (#157) +- `gen_dihedral` keyword argument in `create_residue`, `create_molecule` and `create_hydrogel` to auto-generate dihedral potentials from the bond topology. (#157) +- Dihedral potentials follow the IUPAC/GROMACS convention at the API level: users pass the parameters of the GROMACS proper dihedral (type 1), `V(phi) = k*[1 + cos(n*phi - phi_s)]` with `phi = 0` at the cis conformation, as `k`, `phi_s` and `mult`. The conversion to the ESPResSo `Dihedral` interaction is applied internally by the ESPResSo engine adapter. (#157) +- Sample script `samples/setup_dihedral_potential.py` and test suite `testsuite/dihedral_tests.py` for the new functionality, including tests that compare the energies and forces set up in ESPResSo with the analytical GROMACS expression. (#157) - Simulation-engine abstraction layer to decouple pyMBE molecular-building logic from simulation-engine-specific implementation details. (#151) - Added `EspressoEngine` to delegate ESPResSo-specific operations, including reaction methods, setup utilities, and particle/bond/angle insertion. (#151) - Added `add_instances_to_engine()` to transfer pyMBE particles, bonds, and angles to the selected simulation engine. (#151)