Repository: marrink-lab/vermouth-martinize

Issue:
Incorrect Martini 2.2 topologies for alpha-helical peptides
#### Short description
`vermouth-martinize` generates incorrect Martini 2.2 topologies for alpha-helical peptides. 

#### Detailed description
Using MODELLER, I generated a short 20 amino acids long polyleucine peptide in an alpha-helical conformation. I then used two different versions of `vermouth-martinize` (v0.9.6 and v0.9.7-dev9) and an old version of `martinize.py` script (labeled as version 2.6_3) to coarse-grain the peptide.

Commands used: 
- for `vermouth-martinize`: `-f l20.pdb -o topol.top -x l20_martini.pdb -ss HHHHHHHHHHHHHHHHHHHH -p backbone -ff martini22`,
- for `martinize.py`: `-f l20.pdb -o topol.top -x l20_martini.pdb -ss HHHHHHHHHHHHHHHHHHHH -p backbone`.

Both versions of `vermouth-martinize` generated identical topologies with the following `[ atoms ]` section:

```
[ atoms ]
 1 N0   1 LEU BB   1   1
 2 AC1  1 LEU SC1  2 0.0
 3 N0   2 LEU BB   3 0.0
 4 AC1  2 LEU SC1  4 0.0
 5 N0   3 LEU BB   5 0.0
 6 AC1  3 LEU SC1  6 0.0
 7 N0   4 LEU BB   7 0.0
 8 AC1  4 LEU SC1  8 0.0
 9 N0   5 LEU BB   9 0.0
10 AC1  5 LEU SC1 10 0.0
11 N0   6 LEU BB  11 0.0
12 AC1  6 LEU SC1 12 0.0
13 N0   7 LEU BB  13 0.0
14 AC1  7 LEU SC1 14 0.0
15 N0   8 LEU BB  15 0.0
16 AC1  8 LEU SC1 16 0.0
17 N0   9 LEU BB  17 0.0
18 AC1  9 LEU SC1 18 0.0
19 N0  10 LEU BB  19 0.0
20 AC1 10 LEU SC1 20 0.0
21 N0  11 LEU BB  21 0.0
22 AC1 11 LEU SC1 22 0.0
23 N0  12 LEU BB  23 0.0
24 AC1 12 LEU SC1 24 0.0
25 N0  13 LEU BB  25 0.0
26 AC1 13 LEU SC1 26 0.0
27 N0  14 LEU BB  27 0.0
28 AC1 14 LEU SC1 28 0.0
29 N0  15 LEU BB  29 0.0
30 AC1 15 LEU SC1 30 0.0
31 N0  16 LEU BB  31 0.0
32 AC1 16 LEU SC1 32 0.0
33 N0  17 LEU BB  33 0.0
34 AC1 17 LEU SC1 34 0.0
35 N0  18 LEU BB  35 0.0
36 AC1 18 LEU SC1 36 0.0
37 N0  19 LEU BB  37 0.0
38 AC1 19 LEU SC1 38 0.0
39 N0  20 LEU BB  39  -1
40 AC1 20 LEU SC1 40 0.0
```

`martinize.py` script generated this `[ atoms ]` section:

```
[ atoms ]
    1    Qd     1   LEU    BB     1  1.0000 ; 1
    2   AC1     1   LEU   SC1     2  0.0000 ; 1
    3    Nd     2   LEU    BB     3  0.0000 ; 1
    4   AC1     2   LEU   SC1     4  0.0000 ; 1
    5    Nd     3   LEU    BB     5  0.0000 ; 1
    6   AC1     3   LEU   SC1     6  0.0000 ; 1
    7    Nd     4   LEU    BB     7  0.0000 ; 1
    8   AC1     4   LEU   SC1     8  0.0000 ; 1
    9    N0     5   LEU    BB     9  0.0000 ; H
   10   AC1     5   LEU   SC1    10  0.0000 ; H
   11    N0     6   LEU    BB    11  0.0000 ; H
   12   AC1     6   LEU   SC1    12  0.0000 ; H
   13    N0     7   LEU    BB    13  0.0000 ; H
   14   AC1     7   LEU   SC1    14  0.0000 ; H
   15    N0     8   LEU    BB    15  0.0000 ; H
   16   AC1     8   LEU   SC1    16  0.0000 ; H
   17    N0     9   LEU    BB    17  0.0000 ; H
   18   AC1     9   LEU   SC1    18  0.0000 ; H
   19    N0    10   LEU    BB    19  0.0000 ; H
   20   AC1    10   LEU   SC1    20  0.0000 ; H
   21    N0    11   LEU    BB    21  0.0000 ; H
   22   AC1    11   LEU   SC1    22  0.0000 ; H
   23    N0    12   LEU    BB    23  0.0000 ; H
   24   AC1    12   LEU   SC1    24  0.0000 ; H
   25    N0    13   LEU    BB    25  0.0000 ; H
   26   AC1    13   LEU   SC1    26  0.0000 ; H
   27    N0    14   LEU    BB    27  0.0000 ; H
   28   AC1    14   LEU   SC1    28  0.0000 ; H
   29    N0    15   LEU    BB    29  0.0000 ; H
   30   AC1    15   LEU   SC1    30  0.0000 ; H
   31    N0    16   LEU    BB    31  0.0000 ; H
   32   AC1    16   LEU   SC1    32  0.0000 ; H
   33    Na    17   LEU    BB    33  0.0000 ; 2
   34   AC1    17   LEU   SC1    34  0.0000 ; 2
   35    Na    18   LEU    BB    35  0.0000 ; 2
   36   AC1    18   LEU   SC1    36  0.0000 ; 2
   37    Na    19   LEU    BB    37  0.0000 ; 2
   38   AC1    19   LEU   SC1    38  0.0000 ; 2
   39    Qa    20   LEU    BB    39 -1.0000 ; 2
   40   AC1    20   LEU   SC1    40  0.0000 ; 2
```

Clearly, the bead types for the backbone of the aminoacids near the termini of the peptide differ, with `martinize.py` generating correct bead types. Note especially that  `1 N0   1 LEU BB   1   1` which is a charged backbone bead obviously should not be of type `N0`.

The parameters for bonds, angles, position restraints, and dihedrals seem to be consistent between the individual scripts.

I have attempted to force `vermouth-martinize` to generate the correct topology by providing a slightly different secondary structure. Namely, I used `-ss 1111HHHHHHHHHHHH2222` instead of `-ss HHHHHHHHHHHHHHHHHHHH` but obtained the same (incorrect) result as before.

#### Files
All input and output files as well as the specific version of the `martinize.py` script used can be temporarily found at [github.com/Ladme/vermouth-martinize-bug](https://github.com/Ladme/vermouth-martinize-bug) until the issue is resolved.

#### Acknowledgments
Credit goes to Denys Biriukov and Peter Pajtinka for discovering this issue.

Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index c3fe7be..fae21f5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,6 +2,7 @@ test.pdb
 **/__pycache__
 **.pyc
 **.bak
+\#*
 
 **.lprof
 
diff --git a/doc/source/tutorials/6_adding_residues_links/index.rst b/doc/source/tutorials/6_adding_residues_links/index.rst
index effb5b1..87bd946 100644
--- a/doc/source/tutorials/6_adding_residues_links/index.rst
+++ b/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.
diff --git a/vermouth/data/force_fields/elnedyn21/aminoacids.ff b/vermouth/data/force_fields/elnedyn21/aminoacids.ff
index fd25a19..20adc72 100644
--- a/vermouth/data/force_fields/elnedyn21/aminoacids.ff
+++ b/vermouth/data/force_fields/elnedyn21/aminoacids.ff
@@ -552,22 +552,22 @@ 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
 
@@ -575,19 +575,19 @@ collagen
 [ 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 ]
diff --git a/vermouth/data/force_fields/elnedyn22/aminoacids.ff b/vermouth/data/force_fields/elnedyn22/aminoacids.ff
index b13d4a9..52459be 100644
--- a/vermouth/data/force_fields/elnedyn22/aminoacids.ff
+++ b/vermouth/data/force_fields/elnedyn22/aminoacids.ff
@@ -552,22 +552,22 @@ 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
 
@@ -575,19 +575,19 @@ collagen
 [ 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 ]
diff --git a/vermouth/data/force_fields/elnedyn22p/aminoacids.ff b/vermouth/data/force_fields/elnedyn22p/aminoacids.ff
index 6dfb9ef..dace4d1 100644
--- a/vermouth/data/force_fields/elnedyn22p/aminoacids.ff
+++ b/vermouth/data/force_fields/elnedyn22p/aminoacids.ff
@@ -628,22 +628,22 @@ 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
 
@@ -651,19 +651,19 @@ collagen
 [ 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 ]
diff --git a/vermouth/data/force_fields/martini22/aminoacids.ff b/vermouth/data/force_fields/martini22/aminoacids.ff
index 6d4fe9c..76f7965 100644
--- a/vermouth/data/force_fields/martini22/aminoacids.ff
+++ b/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
diff --git a/vermouth/data/force_fields/martini22p/aminoacids.ff b/vermouth/data/force_fields/martini22p/aminoacids.ff
index 82bfdaa..dc8d1b1 100644
--- a/vermouth/data/force_fields/martini22p/aminoacids.ff
+++ b/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
diff --git a/vermouth/molecule.py b/vermouth/molecule.py
index ea47a2c..cf2429f 100644
--- a/vermouth/molecule.py
+++ b/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
diff --git a/vermouth/processors/do_links.py b/vermouth/processors/do_links.py
index 0017fbb..4284bea 100644
--- a/vermouth/processors/do_links.py
+++ b/vermouth/processors/do_links.py
@@ -25,7 +25,32 @@ from .processor import Processor
 
 
 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):


Test patch:
diff --git a/vermouth/tests/test_links.py b/vermouth/tests/test_links.py
index da2db25..4bd8453 100644
--- a/vermouth/tests/test_links.py
+++ b/vermouth/tests/test_links.py
@@ -18,7 +18,8 @@ import numpy as np
 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