diff --git a/CHANGELOG.md b/CHANGELOG.md
index a15004b..1ecbb02 100644
--- a/CHANGELOG.md
+++ b/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)
diff --git a/pyMBE/pyMBE.py b/pyMBE/pyMBE.py
index 5b0f8ca..d8d7f85 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,248 @@ 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_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:
+ ('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_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"]
+ 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_s": PintQuantity.from_quantity(q=term["phi_s"],
+ 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", 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
+ 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_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:
+ 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)
+ 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_s, 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..d60243d 100644
--- a/pyMBE/simulation_builder/espresso_engine.py
+++ b/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)
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..bf7a0ef
--- /dev/null
+++ b/pyMBE/storage/templates/dihedral.py
@@ -0,0 +1,143 @@
+#
+# 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 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.
+
+ 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_s`` (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")``.
+ - 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")
+ 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..4ae389d
--- /dev/null
+++ b/samples/setup_dihedral_potential.py
@@ -0,0 +1,182 @@
+#
+# 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
+# - 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
+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` (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":
+ 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_s': 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(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")
+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_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,
+ 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[1]))} ESPResSo Dihedral objects "
+ f"for a {len(multi_term_parameters)}-term dihedral")
+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:
+ 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) 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_s = np.deg2rad(gromacs_phi_s_deg)
+
+pmb.define_dihedral_potential(dihedral_type="periodic",
+ dihedral_parameters={'k': gromacs_k * pmb.units('reduced_energy'),
+ 'phi_s': phi_s * 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[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_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_s + np.pi))
+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..68854f0 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_s": 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..832aa14
--- /dev/null
+++ b/testsuite/dihedral_tests.py
@@ -0,0 +1,933 @@
+#
+# 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
+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):
+ 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_s': np.pi * pmb.units(''),
+ 'mult': 1,
+ }
+
+ self.periodic_dihedral_params_alt = {
+ 'k': 3 * pmb.units('reduced_energy'),
+ 'phi_s': np.pi / 2 * pmb.units(''),
+ 'mult': 2,
+ }
+
+ self.multi_term_dihedral_params = [
+ {'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 = {
+ 'r_0': 0.4 * pmb.units.nm,
+ 'k': 400 * pmb.units('reduced_energy / reduced_length**2'),
+ }
+
+ def get_dihedral_objects(self, particle_id2):
+ """
+ Returns the list of dihedral interaction objects stored in ESPResSo for a
+ 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_id2).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.
+ pyMBE takes GROMACS type 1 parameters; the ESPResSo phase is -(phi_s + pi).
+
+ 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_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)
+ 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_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])
+
+ 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_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)
+ 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_id2=pids[1])
+ 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_id2=pids[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_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_s': 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_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_s': 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_id2=pids[1])
+ 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_id2=pids[1])), 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")
+
+ 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/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()
diff --git a/testsuite/test_io_database.py b/testsuite/test_io_database.py
index 7667f36..1768187 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_s": 1.0 * pmb.units(""),
+ "mult": 1}
+ 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"),
+ ("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_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]
+ 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