Merge pull request #561 from marrink-lab/issue_560

Issue 560

.gitignore

@@ -2,6 +2,7 @@ test.pdb
 **/__pycache__
 **.pyc
 **.bak
+\#*
 
 **.lprof
 

doc/source/tutorials/6_adding_residues_links/index.rst

@@ -270,3 +270,18 @@ Which will produce the exact same topology.
 If you *do* need to add a residue that can be used in any kind of protein
 please take a look at how the Martini 3 force field is implemented, and deals
 with e.g. the secondary structure dependence.
+
+Links and Modifications
++++++++++++++++++++++++
+Something to keep in mind is that Links get applied after Modifications (at the time of writing). This can mean that
+your Link overwrites, for example, terminal parameters. For this reason, you can filter nodes where Links get applied
+much like you can limit Links by atom names or secondary structure. In particular, you can add a ``"modifications"``
+attribute to links nodes. This follows the following rules:
+
+1. Links that don't specify modifications simply match.
+2. Links that specify empty modifications (``"modifications": []`` or ``null``) only match atoms that have no
+   modifications.
+3. Links that specify a list of modifications (``"modifications": ["C-ter", "ASP-HD2"]``) only match atoms that carry
+   that exact set of modifications.
+4. Links that specify a string or Choice of modifications (``"modifications": "C-ter"`` or ``"C-ter|COOH-ter"``) only
+   match atoms where *all* the atoms modifications match.

vermouth/data/force_fields/elnedyn21/aminoacids.ff

@@ -552,42 +552,42 @@ TRP                1
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}}
+BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}}
+BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}, "modifications": null}
 [ features ]
 collagen
 
 ;; Fix bead types for ALA and PRO.
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}, "modifications": []}
 
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}, "modifications": []}
 [ features ]
 collagen
 
 [ link ]
 resname "PRO"
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": []}
 
 ;; Setup the bonds.
 [ link ]

vermouth/data/force_fields/elnedyn22/aminoacids.ff

@@ -552,42 +552,42 @@ TRP                1
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}}
+BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}}
+BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}, "modifications": null}
 [ features ]
 collagen
 
 ;; Fix bead types for ALA and PRO.
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}, "modifications": null}
 
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}, "modifications": null}
 [ features ]
 collagen
 
 [ link ]
 resname "PRO"
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 ;; Setup the bonds.
 [ link ]

vermouth/data/force_fields/elnedyn22p/aminoacids.ff

@@ -628,42 +628,42 @@ TRP                1
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}}
+BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}}
+BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}, "modifications": null}
 [ features ]
 collagen
 
 ;; Fix bead types for ALA and PRO.
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}, "modifications": null}
 
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}, "modifications": null}
 [ features ]
 collagen
 
 [ link ]
 resname "PRO"
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 ;; Setup the bonds.
 [ link ]

vermouth/data/force_fields/martini22/aminoacids.ff

@@ -621,45 +621,45 @@ BB -BB
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}}
+BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}}
+BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ features ]
 collagen
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}, "modifications": null}
 
 ;; Fix bead types for ALA and PRO.
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "S", "replace": {"atype": "P4"}}
+BB {"cgsecstruct": "S", "replace": {"atype": "P4"}, "modifications": null}
 
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}, "modifications": null}
 
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}, "modifications": null}
 
 [ link ]
 resname "PRO"
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 ;; Setup the bonds. We only have the bonds assuming everything is coil.
 ;; We always select the lowest force constant when the two residues involved

vermouth/data/force_fields/martini22p/aminoacids.ff

@@ -644,43 +644,43 @@ BB -BB
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}}
+BB {"cgsecstruct": "T|3|E", "replace": {"atype": "Nda"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}}
+BB {"cgsecstruct": "1", "replace": {"atype": "Nd"}, "modifications": null}
 
 [ link ]
 resname $protein_resnames
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "N0"}, "modifications": null}
 
 ;; Fix bead types for ALA and PRO.
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "S", "replace": {"atype": "P4"}}
+BB {"cgsecstruct": "S", "replace": {"atype": "P4"}, "modifications": null}
 
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}}
+BB {"cgsecstruct": "T|3|2|1|E", "replace": {"atype": "N0"}, "modifications": null}
 
 [ link ]
 resname "ALA|PRO|HYP"
 [ atoms ]
-BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}}
+BB {"cgsecstruct": "H|F", "replace": {"atype": "C5"}, "modifications": null}
 
 [ link ]
 resname "PRO"
 [ atoms ]
-BB {"cgsecstruct": "2", "replace": {"atype": "Na"}}
+BB {"cgsecstruct": "2", "replace": {"atype": "Na"}, "modifications": null}
 
 ;; Setup the bonds. We only have the bonds assuming everything is coil.
 ;; We always select the lowest force constant when the two residues involved

vermouth/molecule.py

@@ -1350,7 +1350,7 @@ def attributes_match(attributes, template_attributes, ignore_keys=()):
 
     Returns ``True`` if the attributes from the link match the ones from the
     molecule; returns ``False`` otherwise. The attributes from a link match
-    with those of a molecule is all the individual attribute from the link
+    with those of a molecule if all the individual attribute from the link
     match the corresponding ones in the molecule. In the simplest case, these
     attribute match if their values are equal. If the value of the link
     attribute is an instance of :class:`LinkPredicate`, then the attributes

vermouth/processors/do_links.py

@@ -25,7 +25,32 @@
 
 
 def _atoms_match(node1, node2):
-    return attributes_match(node1, node2, ignore_keys=('order', 'replace'))
+    # node1 is molecule, node2 is link
+    # modifications are named as tuples, e.g. `('C-ter',)`, since mappings can
+    # deal with multiple modifications at the same time. Here we build a flat
+    # list of relevant modification names, and require that /all/ of these match
+    # the link['modifications']
+    mods = []
+    for mod in node1.get('modifications', []):
+        mods.extend(mod.name)
+
+                  # No modifications specified by link: always match
+    mods_match = ('modifications' not in node2 or
+                  # empty modifications in link and no modifications in molecule: match
+                  (not node2['modifications'] and not mods) or
+                  # Else, if both specify modifications, then...
+                  (node2['modifications'] and mods and
+                   # link modifications must be a list, and molecule mods must
+                   # match links mods exactly
+                   ((isinstance(node2['modifications'], list) and sorted(mods) == sorted(node2['modifications'])) or
+                    # Or link mods are a simple string or a Choice, and all
+                    # molecule modifications must be accounted for
+                    # Here we need to do a little jiggery-pokery to leverage
+                    # attributes_match. This probably means that that function
+                    # needs to be cut up into smaller pieces.
+                    all(attributes_match({'_': modname}, {'_': node2['modifications']}) for modname in mods))))
+
+    return bool(mods_match and attributes_match(node1, node2, ignore_keys=('order', 'replace', 'modifications')))
 
 
 def _is_valid_non_edges(molecule, link, rev_raw_match):

vermouth/tests/test_links.py

@@ -18,7 +18,8 @@
 import pytest
 
 from vermouth.processors import do_links, DoLinks
-from vermouth.molecule import Molecule, Link
+from vermouth.molecule import Molecule, Link, Modification, Choice
+from vermouth.processors.do_links import _atoms_match
 import vermouth.forcefield
 
 
@@ -278,3 +279,38 @@ def test_link_processor(mol_nodes, mol_edges, link_nodes, link_edges,
     out = DoLinks().run_molecule(mol)
     assert dict(out.nodes(data=True)) == dict(expected_nodes)
     assert set(out.edges(data=False)) == set(expected_edges)
+
+
+def modification(*names):
+    mod = Modification(name=tuple(names))
+    return mod
+
+
+@pytest.mark.parametrize('mol_node, link_node, expected', (
+    [{}, {}, True],  # No restrictions, match
+    [{'modifications': [modification('A')]}, {}, True],  # No restrictions, match
+    [{}, {'modifications': []}, True],  # No link modifications, no mol modifications, match
+    [{}, {'modifications': None}, True],  # No link modifications, no mol modifications, match
+    [{'modifications': [modification('A')]}, {'modifications': None}, False],  # No mods allowed, no match
+    [{'modifications': [modification('A')]}, {'modifications': []}, False],  # No mods allowed, no match
+    [{'modifications': []}, {'modifications': []}, True],  # No link modifications, no mol modifications, match
+    [{'modifications': []}, {'modifications': None}, True],  # No link modifications, no mol modifications, match
+    [{}, {'modifications': ['A']}, False],   # Link modifications, but no mol modifications, don't match
+    [{'modifications': []}, {'modifications': 'A'}, False],  # Link modifications, but no mol modifications, don't match
+    [{'modifications': [modification('A')]}, {'modifications': 'A'}, True],  # Link modifications and matching mol mods, match
+    [{'modifications': [modification('A')]}, {'modifications': 'B'}, False],  # Link modifications but no matching mol mods, don't match
+    [{'modifications': [modification('A', 'B')]}, {'modifications': 'B'}, False],  # Link modifications but no not all mol mods match, don't match
+    [{'modifications': [modification('A'), modification('B')]}, {'modifications': 'B'}, False],  # Link modifications but no not all mol mods match, don't match
+    [{'modifications': [modification('A')]}, {'modifications': Choice(['A', 'B'])}, True],  # Link modifications and matching mol mods, match
+    [{'modifications': [modification('A', 'B')]}, {'modifications': ['A', 'B']}, True],  # Link modifications and matching mol mods, match
+    [{'modifications': [modification('A', 'B')]}, {'modifications': ['B', 'A']}, True],  # Link modifications and matching mol mods, match
+    [{'modifications': [modification('A'), modification('B')]}, {'modifications': ['A', 'B']}, True],  # Link modifications and matching mol mods, match
+    [{'modifications': [modification('A'), modification('B')]}, {'modifications': 'A'}, False],  # Unmatched mod, no match
+    [{'modifications': [modification('A'), modification('A')]}, {'modifications': 'A'}, True],  # Multiple mods, all matched
+    [{'modifications': [modification('A'), modification('B')]}, {'modifications': Choice(['A', 'B'])}, True],  # Choice
+    [{'modifications': [modification('A'), modification('A')]}, {'modifications': Choice(['A', 'B'])}, True],  # Choice
+    [{'modifications': [modification('C'), modification('A')]}, {'modifications': Choice(['A', 'B'])}, False],  # Unmatched choice
+))
+def test_modification_matching(mol_node, link_node, expected):
+    found = _atoms_match(mol_node, link_node)
+    assert found == expected