Add support for dihedral potentials

CHANGELOG.md

@@ -8,6 +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. (#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)

pyMBE/simulation_builder/espresso_engine.py

@@ -57,7 +57,30 @@ 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)
+
+        # 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_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',
+                                    value=True)
+
     def _add_bond(self,particle_id1,particle_id2,bond_inst):
         """
         Helper function to add bond instances to ESPResSo.
@@ -248,7 +271,70 @@ 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_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
+
+        if dihedral_type == "periodic":
+            dihedral_insts = []
+            for term in dihedral_terms:
+                k = term["k"].m_as("reduced_energy")
+                phi_s = float(term["phi_s"].magnitude)
+                mult = int(round(float(term["mult"].magnitude)))
+                # 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")
+        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 +1440,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 +1499,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)
 

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 <http://www.gnu.org/licenses/>.
+#
+
+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

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,