Repository: marrink-lab/vermouth-martinize
Issue: Selectively removing interactions from molecule

### Discussed in https://github.com/marrink-lab/vermouth-martinize/discussions/482

Regarding reading in itps: it would be useful to have read back in header information and meta attributes for different bond groups. The former would be useful for checking aspects of secondary structure assignment, and the latter for itp editing when needed.

<div type='discussions-op-text'>

<sup>Originally posted by **csbrasnett** November 16, 2022</sup>
I need to edit an itp to remove parts of the elastic network, as determined by secondary structure assignment.

As a minimal example, I would have thought something like this should work:

```python
# loop over regions where we want to keep EN
for i in wanted_region:
     # look at all the elastic network interactions
     for j in all_EN_interactions:
         # compare atom indicies of the EN interaction to the boundaries of the region we want to keep
         # if they're within, keep them
         if all_EN_interactions[j] in wanted_region[i]:
             continue
         # if they're outside, remove the interaction
         else:
             ff.blocks[molname].remove_interaction('bonds', atoms = all_EN_interactions[j])
```

However, I keep getting a KeyError
`KeyError: "Can't find interaction of type bonds between atoms ['0', '5'] and with version 0"`
even though I can definitely see that there's an interaction there. I've tried various versions of the atoms argument in remove_interaction (give as an array, a list, a list of strings, etc.) but none of these seem to work.

I guess it might want me to give a version explicitly, but this doesn't seem to work.

Perhaps this should be raised as a separate issue, but on this last point when I read in the itp file that's been generated with martinize2, the meta information isn't read back in by vermouth - I think this might be quite a useful feature to have for this sort of purpose?



</div>

Changed files:
vermouth/molecule.py, vermouth/tests/test_molecule.py

diff --git a/vermouth/molecule.py b/vermouth/molecule.py
index b4f9bd0..af1d41c 100644
--- a/vermouth/molecule.py
+++ b/vermouth/molecule.py
@@ -589,13 +589,15 @@ class Molecule(nx.Graph):
         """
         idx = 0
         for idx, interaction in enumerate(self.interactions[type_]):
-            if interaction.atoms == atoms and interaction.meta.get('version', 0):
+            if interaction.atoms == atoms and interaction.meta.get('version', 0) == version:
                 break
         else:  # no break
             msg = ("Can't find interaction of type {} between atoms {} "
                    "and with version {}")
             raise KeyError(msg.format(type_, atoms, version))
         del self.interactions[type_][idx]
+        if not self.interactions[type_]:
+            del self.interactions[type_]
 
     def remove_matching_interaction(self, type_, template_interaction):
         """
diff --git a/vermouth/tests/test_molecule.py b/vermouth/tests/test_molecule.py
index 57374fc..e395431 100644
--- a/vermouth/tests/test_molecule.py
+++ b/vermouth/tests/test_molecule.py
@@ -1279,3 +1279,58 @@ def test_to_molecule():
                              parameters=['a', '0.1', '300'],
                              meta={'b': 1}),]
     assert ref_bonds == test_molecule.interactions['bonds']
+
+
+@pytest.mark.parametrize('atoms, bonds, interactions, removed, expected', [
+    # empty molecule
+    ([], [], [], {}, {}),
+    # interactions that all need to be removed
+    ([1, 2], [(1, 2)], [('bond', (1, 2), {})], {'type_': 'bond', "atoms": (1, 2)}, {}),
+    # Molecule with interactions of different types of which some need to be removed
+    (
+        [1, 2, 3, 4, 5, 6, 7],
+        [(1, 2), (3, 4), (5, 7), (6, 7)],
+        [('bond', (1, 6), {}), ('bond', (5, 4), {}), ('angle', (1, 2), {})],
+        {'type_': 'bond', 'atoms': (5, 4)},
+        {'angle': [vermouth.molecule.Interaction(atoms=(1, 2), meta={}, parameters={})],
+         'bond': [vermouth.molecule.Interaction(atoms=(1, 6), meta={}, parameters={})]}
+    ),
+    (
+        [1, 2, 3, 4],
+        [(1, 2), (2, 3), (3, 4)],
+        [('bond', (1, 2), {}), ('bond', (2, 3), {}), ('bond', (2, 3), {}, {'version': 1})],
+        {'type_': 'bond', 'atoms': (2, 3), 'version': 1},
+        {'bond': [
+            vermouth.molecule.Interaction(atoms=(1, 2), meta={}, parameters={}),
+            vermouth.molecule.Interaction(atoms=(2, 3), meta={}, parameters={}),
+        ]},
+    ),
+    (
+            [1, 2, 3, 4, 5, 6, 7],
+            [(1, 2), (3, 4), (5, 7), (6, 7)],
+            [('bond', (1, 6), [1, 2, 3], {'comment': 'a'}), ('bond', (5, 4), {'comment': 'b'}), ('angle', (1, 2), (4, 5, 6), {'version': 1})],
+            {'type_': 'bond', 'atoms': (5, 4)},
+            {'angle': [vermouth.molecule.Interaction(atoms=(1, 2), meta={'version': 1}, parameters=(4, 5, 6))],
+             'bond': [vermouth.molecule.Interaction(atoms=(1, 6), meta={'comment': 'a'}, parameters=[1, 2, 3])]}
+    ),
+
+])
+def test_remove_interaction(atoms, bonds, interactions, removed, expected):
+    """
+    Test whether molecule.remove_node also removes the corresponding
+    interactions
+    """
+    molecule = vermouth.molecule.Molecule()
+    molecule.add_nodes_from(atoms)
+    molecule.add_edges_from(bonds)
+    for type_, atoms, params, *meta in interactions:
+        if meta:
+            meta = meta[0]
+        else:
+            meta = None
+        molecule.add_interaction(type_, atoms, params, meta=meta)
+
+    if removed:
+        molecule.remove_interaction(**removed)
+
+    assert molecule.interactions == expected
\\ No newline at end of file