Merge pull request #166 from marrink-lab/fix_bug_warning

skip unused molecules in residue check
marrink-lab · Nov 16, 2021 · e95b691 · e95b691
2 parents 71910f2 + 7edad65
commit e95b691
Show file tree

Hide file tree

Showing 2 changed files with 66 additions and 0 deletions.
diff --git a/polyply/src/check_residue_equivalence.py b/polyply/src/check_residue_equivalence.py
@@ -20,6 +20,10 @@ def check_residue_equivalence(topology):
     resids = {}
     molnames = {}
     for mol_name, mol_idxs in topology.mol_idx_by_name.items():
+        # molecules that are included but not used in the topology
+        # simply don't get mol_idxs but we need to skip them here
+        if not mol_idxs:
+            continue
         molecule = topology.molecules[mol_idxs[0]]
         for node in molecule.nodes:
             resname = molecule.nodes[node]["resname"]

diff --git a/polyply/tests/test_residue_check.py b/polyply/tests/test_residue_check.py
@@ -128,3 +128,65 @@ def test_raise_residue_error(top_lines):
     topology.volumes = {"GLY": 0.53, "GLU": 0.67, "ASP": 0.43}
     with pytest.raises(IOError):
         polyply.src.check_residue_equivalence.check_residue_equivalence(topology)
+
+def test_raise_residue_no_error():
+    """
+    This test makes sure that we actually skip molecules that are defined
+    in the top file but not used in the actual system.
+    """
+    top_lines="""
+    [ defaults ]
+    1   1   no   1.0     1.0
+    [ atomtypes ]
+    N0 72.0 0.000 A 0.0 0.0
+    N1 72.0 0.000 A 0.0 0.0
+    [ nonbond_params ]
+    N0   N0   1 4.700000e-01    3.700000e+00
+    N1   N1   1 4.700000e-01    3.700000e+00
+    [ moleculetype ]
+    testA 1
+    [ atoms ]
+    1    N0   1   GLY    BB   1 0.00     45
+    2    N0   1   GLY    SC1  1 0.00     45
+    3    N0   1   GLY    SC2  1 0.00     45
+    4    N1   2   GLU    BB   2 0.00     45
+    5    N1   2   GLU    SC1  2 0.00     45
+    6    N1   2   GLU    SC2  2 0.00     45
+    7    N1   2   GLU    SC3  2 0.00     45
+    [ bonds ]
+    1     2    1  0.47 2000
+    2     3    1  0.47 2000
+    3     4    1  0.47 2000
+    4     5    1  0.47 2000
+    5     6    1  0.47 2000
+    6     7    1  0.47 2000
+    [ moleculetype ]
+    testB 1
+    [ atoms ]
+    1    N0   1   ASP    BB   1 0.00     45
+    2    N0   1   ASP    SC3  1 0.00     45
+    3    N0   1   ASP    SC4  1 0.00     45
+    4    N1   2   GLU    BB   2 0.00     45
+    5    N1   2   GLU    SC1  2 0.00     45
+    6    N1   2   GLU    SC2  2 0.00     45
+    7    N1   2   GLU    SC3  2 0.00     45
+    [ bonds ]
+    1     2    1  0.47 2000
+    2     3    1  0.47 2000
+    3     4    1  0.47 2000
+    4     5    1  0.47 2000
+    5     6    1  0.47 2000
+    5     7    1  0.47 2000
+    [ system ]
+    test system
+    [ molecules ]
+    testA 1
+    """
+    lines = textwrap.dedent(top_lines)
+    lines = lines.splitlines()
+    force_field = ForceField("test")
+    topology = Topology(force_field)
+    read_topology(lines=lines, topology=topology, cwdir="./")
+    topology.preprocess()
+    topology.volumes = {"GLY": 0.53, "GLU": 0.67, "ASP": 0.43}
+    polyply.src.check_residue_equivalence.check_residue_equivalence(topology)