Merge top

polyply/src/topology.py

@@ -28,6 +28,9 @@
 from .top_parser import read_topology
 from .linalg_functions import center_of_geometry
 
+class MergeError(Exception):
+    """ Raised when two objects cannot be merged"""
+
 COORD_PARSERS = {"pdb": read_pdb,
                  "gro": read_gro}
 
@@ -192,7 +195,7 @@ def match_dihedral_interaction_types(atoms, interaction_dict):
 class Topology(System):
     """
     Ties together vermouth molecule definitions, and
-    Gromacs topology information.
+    Gromacs or other programs topology information.
 
     Parameters
     ----------
@@ -212,21 +215,86 @@ class Topology(System):
         A dictionary of all typed parameter
     defines: list
         A list of everything that is defined
+    defaults: dict
+        A dict with default entries of GMX topfile
+    atom_types: dict
     """
 
     def __init__(self, force_field, name=None):
+        # self.molecules with super
+        # self.foce_field with super
         super().__init__(force_field)
         self.name = name
+        # simple dict attributes
         self.defaults = {}
         self.defines = {}
-        self.description = []
         self.atom_types = {}
-        self.types = defaultdict(dict)
         self.nonbond_params = {}
-        self.mol_idx_by_name = defaultdict(list)
+        self.volumes = {}
+        # list attributes
+        self.description = []
         self.persistences = []
+        self.molecules = []
+        # default dict list
+        self.mol_idx_by_name = defaultdict(list)
+        # default_dict dict
         self.distance_restraints = defaultdict(dict)
-        self.volumes = {}
+        self.types = defaultdict(dict)
+
+    def merge(self, other_top, check_duplicates=True):
+        """
+        Merge two topologies updating their attribute dictionaries.
+        If `check_duplicates` is True attribute value combinations
+        that already exist will be overwritten by those in
+        the to be merged topology. An error will raise otherwise if
+        there are any conflicting duplicate entries.
+
+        Note distance restraints are not checked for possible
+        conflicting definitions.
+
+        Parameters
+        ----------
+        other_top: :class:`polyply.src.topology.Topology`
+            the topology to be merged
+
+        check_duplicates: bool
+            check for conflicting attribute value pairs
+
+        """
+        for attribute in ["defaults", "defines", "atom_types", "nonbond_params", "volumes"]:
+            self_dict = getattr(self, attribute)
+            if check_duplicates:
+                for key, value in getattr(other_top, attribute).items():
+                    if key in self_dict and self_dict[key] != value:
+                        msg = f"Conflicting entry in {attribute} with key {key}"
+                        raise MergeError(msg)
+            getattr(self, attribute).update(getattr(other_top, attribute))
+
+        self.description += other_top.description
+        self.persistences += other_top.persistences
+
+        for atoms, dist_restr in other_top.distance_restraints:
+            self.distance_restraints[atoms].append(dist_restr)
+
+        for molecule in other_top.molecules:
+            self.append_molecules(molecule, molecule.mol_name)
+
+        for inter_type in other_top.types:
+            for atoms, type_params in other_top.types[inter_type].items():
+                if atoms in self.types[inter_type] and self.types[inter_type][atoms] != type_params:
+                    typestring = inter_type[:-1]
+                    msg = f"Conflicting entry in {typestring}types for atoms {atoms}"
+                    raise MergeError(msg)
+                self.types[inter_type][atoms] = type_params
+
+    def append_molecules(self, new_molecule, molname):
+        """
+        Add a new molecule to molecule list and update
+        the mol_idx_by_name dictionary.
+        """
+        mol_idx = len(self.molecules)
+        self.molecules.append(new_molecule)
+        self.mol_idx_by_name[molname].append(mol_idx)
 
     def preprocess(self):
         """
@@ -429,6 +497,11 @@ def add_positions_from_file(self, path, skip_res=[], resolution='mol'):
                     meta_mol.nodes[meta_node]["backmap"] = False
 
     def convert_to_vermouth_system(self):
+        """
+        Create a vermouth system by setting the force-field
+        attribute and the molecule list appending the molecule
+        part of the meta-molecules.
+        """
         system = System()
         system.molecules = []
         system.force_field = self.force_field

polyply/tests/test_topology.py

@@ -25,6 +25,12 @@
 from polyply import TEST_DATA
 from polyply.src.topology import Topology
 
+def _lines_to_topology(top, lines):
+    new_lines = textwrap.dedent(lines)
+    new_lines = new_lines.splitlines()
+    polyply.src.top_parser.read_topology(new_lines, top)
+    return top
+
 class TestTopology:
 
     @staticmethod
@@ -652,3 +658,197 @@ def test_replace_types_fail():
         polyply.src.top_parser.read_topology(new_lines, top)
         with pytest.raises(OSError):
              top.gen_bonded_interactions()
+
+    @staticmethod
+    @pytest.mark.parametrize('target_lines, other_lines, method', (
+        (
+        """
+        [ defaults ]
+        1.0   1.0   yes  1.0     1.0
+        [ bondtypes ]
+        C       C       1       0.1335  502080.0
+        [ moleculetype ]
+        test 3
+        [ atoms ]
+        1 C   1 test C1 1   0.0 14.0
+        2 C   1 test C2 2   0.0 12.0
+        [ bonds ]
+        1 2  1
+        [ system ]
+        some title
+        [ molecules ]
+        test 1
+        """,
+        """
+        [ defaults ]
+        1.0   1.0   yes  1.0     1.0
+        [ bondtypes ]
+        P       P       1       0.1335  502080.0
+        [ moleculetype ]
+        other 3
+        [ atoms ]
+        1 P   1 other  C1 1   0.0 14.0
+        2 P   1 other  C2 2   0.0 12.0
+        [ bonds ]
+        1 2  1
+        [ system ]
+        other title
+        [ molecules ]
+        other 1
+        """,
+        True
+        ),
+        ("""
+        [ defaults ]
+        1.0   1.0   yes  1.0     1.0
+        [ bondtypes ]
+        C       C       1       0.1335  502080.0
+        [ moleculetype ]
+        test 3
+        [ atoms ]
+        1 C   1 test C1 1   0.0 14.0
+        2 C   1 test C2 2   0.0 12.0
+        [ bonds ]
+        1 2  1
+        [ system ]
+        some title
+        [ molecules ]
+        test 1
+        """,
+        """
+        [ defaults ]
+        1.0   1.0   yes  1.0     1.0
+        [ bondtypes ]
+        P       P       1       0.1335  502080.0
+        [ moleculetype ]
+        other 3
+        [ atoms ]
+        1 P   1 other  C1 1   0.0 14.0
+        2 P   1 other  C2 2   0.0 12.0
+        [ bonds ]
+        1 2  1
+        [ system ]
+        other title
+        [ molecules ]
+        other 1
+        """,
+        False
+        ),
+        ))
+    def test_merge_topologies(target_lines, other_lines, method):
+        force_field = vermouth.forcefield.ForceField(name='test_ff')
+        target_top =  Topology(force_field, name="target")
+        n_mols_target = len(target_top.molecules)
+        other_top =  Topology(force_field, name="other")
+        _lines_to_topology(target_top, target_lines)
+        _lines_to_topology(other_top, other_lines)
+        target_top.merge(other_top, method)
+        # check if the attributes have both topologies values
+        for attribute in ["defaults", "defines", "atom_types", "nonbond_params"]:
+            target_dict = getattr(target_top, attribute)
+            other_dict = getattr(other_top, attribute)
+            for key in other_dict:
+                assert other_dict[key] == target_dict[key]
+        # check that all descriptions and persistences are in the target
+        # topology
+        for descr in other_top.description:
+            assert descr in target_top.description
+
+        for persis in other_top.persistences:
+            assert persis in target_top.persistences
+
+        for atoms, dist_restr in other_top.distance_restraints:
+            assert dist_restr in self.distance_restraints[atoms]
+
+        for idx, molecule in enumerate(other_top.molecules):
+            assert list(molecule.nodes) == list(target_top.molecules[idx+n_mols_target].nodes)
+            assert list(molecule.edges) == list(target_top.molecules[idx+n_mols_target].edges)
+
+        for inter_type in other_top.types:
+            for atoms, type_params in other_top.types[inter_type].items():
+                assert target_top.types[inter_type][atoms] == type_params
+
+    @staticmethod
+    @pytest.mark.parametrize('target_lines, other_lines', (
+        (
+        """
+        [ defaults ]
+        1.0   1.0   yes  1.0     1.0
+        [ bondtypes ]
+        C       C       1       0.1335  502080.0
+        [ moleculetype ]
+        test 3
+        [ atoms ]
+        1 C   1 test C1 1   0.0 14.0
+        2 C   1 test C2 2   0.0 12.0
+        [ bonds ]
+        1 2  1
+        [ system ]
+        some title
+        [ molecules ]
+        test 1
+        """,
+        """
+        [ defaults ]
+        1.0   2.0   yes  1.0     1.0
+        [ bondtypes ]
+        P       P       1       0.1335  502080.0
+        [ moleculetype ]
+        other 3
+        [ atoms ]
+        1 P   1 other  C1 1   0.0 14.0
+        2 P   1 other  C2 2   0.0 12.0
+        [ bonds ]
+        1 2  1
+        [ system ]
+        other title
+        [ molecules ]
+        other 1
+        """,
+        ),
+        (
+        """
+        [ defaults ]
+        1.0   1.0   yes  1.0     1.0
+        [ bondtypes ]
+        C       C       1       0.1335  502080.0
+        [ moleculetype ]
+        test 3
+        [ atoms ]
+        1 C   1 test C1 1   0.0 14.0
+        2 C   1 test C2 2   0.0 12.0
+        [ bonds ]
+        1 2  1
+        [ system ]
+        some title
+        [ molecules ]
+        test 1
+        """,
+        """
+        [ defaults ]
+        1.0   1.0   yes  1.0     1.0
+        [ bondtypes ]
+        P       P       1       0.1335  502080.0
+        C       C       1       0.13  502080.0
+        [ moleculetype ]
+        other 3
+        [ atoms ]
+        1 P   1 other  C1 1   0.0 14.0
+        2 P   1 other  C2 2   0.0 12.0
+        [ bonds ]
+        1 2  1
+        [ system ]
+        other title
+        [ molecules ]
+        other 1
+        """,
+        ), ))
+    def test_merge_topologies_fail(target_lines, other_lines):
+        force_field = vermouth.forcefield.ForceField(name='test_ff')
+        target_top =  Topology(force_field, name="target")
+        n_mols_target = len(target_top.molecules)
+        other_top =  Topology(force_field, name="other")
+        _lines_to_topology(target_top, target_lines)
+        _lines_to_topology(other_top, other_lines)
+        with pytest.raises(polyply.src.topology.MergeError):
+            target_top.merge(other_top)