raise warning if geometry is out of bounds

polyply/src/build_file_parser.py

@@ -51,10 +51,16 @@ def _molecule(self, line, lineno=0):
         tokens = line.split()
         self.current_molname = tokens[0]
         self.current_molidxs = np.arange(float(tokens[1]), float(tokens[2]), 1., dtype=int)
+        # raise error if molname not in system
+        if tokens[0] not in self.topology.mol_idx_by_name:
+            msg = "Molecule with name {name} is not part of the system."
+            raise IOError(msg.format(name=tokens[0]))
+
         for idx in self.current_molidxs:
             if idx not in self.topology.mol_idx_by_name[tokens[0]]:
-                LOGGER.warning("parsing build file: could not find molecule with name {name} and index {index}.",
-                              **{"index": idx, "name": tokens[0]})
+                msg = ("parsing build file: could not find molecule with "
+                       "name {name} and index {index}.")
+                LOGGER.warning(msg, index=idx, name=tokens[0])
 
     @SectionLineParser.section_parser('molecule', 'cylinder', geom_type="cylinder")
     @SectionLineParser.section_parser('molecule', 'sphere', geom_type="sphere")
@@ -207,7 +213,6 @@ def finalize(self, lineno=0):
         if that molecule is mentioned in the build file by name.
         """
         for mol_idx, molecule in enumerate(self.molecules):
-
             if (molecule.mol_name, mol_idx) in self.build_options:
                 for option in self.build_options[(molecule.mol_name, mol_idx)]:
                     self._tag_nodes(molecule, "restraints", option, molecule.mol_name)
@@ -231,17 +236,15 @@ def _tag_nodes(molecule, keyword, option, molname=""):
                                                 [option['parameters']]
             # broadcast warning if we find the resid but it doesn't match the resname
             elif molecule.nodes[node]["resid"] in resids and not\
-                 molecule.nodes[node]["resname"] == resname:
-                 msg = "parsing build file: could not find resid {resid} with resname {resname} in molecule {molname}."
-                 LOGGER.warning(msg, **{"resid": molecule.nodes[node]["resid"], "resname": resname,
-                                          "molname": molname})
+                molecule.nodes[node]["resname"] == resname:
+                msg = "parsing build file: could not find resid {resid} with resname {resname} in molecule {molname}."
+                LOGGER.warning(msg, resid=molecule.nodes[node]["resid"], resname=resname, molname=molname)
 
             # broadcast warning if we find the resname but it doesn't match the resid
             elif molecule.nodes[node]["resname"] == resname and not\
-                 molecule.nodes[node]["resid"]:
-                 msg = "parsing build file: could not find residue {resname} with resid {resid} in molecule {molname}."
-                 LOGGER.warning(msg, **{"resid": molecule.nodes[node]["resid"], "resname": resname,
-                                        "molname": molname})
+                molecule.nodes[node]["resid"] in resids:
+                msg = "parsing build file: could not find residue {resname} with resids {resids} in molecule {molname}."
+                LOGGER.warning(msg, resids=",".join(map(str, resids)), resname=resname, molname=molname)
 
     @staticmethod
     def _base_parser_geometry(tokens, _type):
@@ -258,8 +261,53 @@ def _base_parser_geometry(tokens, _type):
 
         parameters.append(_type)
         geometry_def["parameters"] = parameters
+
+        is_good_geometry = _check_geometry_def(geometry_def, _type)
+        if not is_good_geometry:
+            msg = ("Geometry restriction {_type} {resname} {start} {stop} "
+                   "extends beyond definite positive coordinates relative "
+                   "to the center point {x:.3f} {y:.3f} {z:.3f}. Be aware "
+                   "polyply only builds in positive coordinate space.")
+            LOGGER.warning(msg, _type=_type, x=point[0], y=point[1], z=point[2],
+                           resname=tokens[0], start=tokens[1], stop=tokens[2])
         return geometry_def
 
+def _check_geometry_def(geom_def, geom_type):
+    """
+    Raise a warning if the point of reference
+    for the geometry type is too close to the
+    origin.
+
+    Parameters
+    ----------
+    geom_def: dict
+        dict with entries: resname, start, stop, point, parameters
+    geom_type: str
+        one of sphere, cylinder, rectangle
+    """
+    msg = ("Geometry restriction {geom_def} extends beyond definite "
+           "positive coordinates relative to the center point "
+           "{x:.3f} {y:.3f} {z:.3f}. Be aware polyply only builds "
+           "in positive coordinate space.")
+
+    point = geom_def['parameters'][1]
+    if geom_type == "sphere":
+        radius = geom_def['parameters'][2]
+        if any( i >= j for i, j in zip([radius, radius, radius], point)):
+            return False
+
+    if geom_type == "rectangle":
+        a, b, c = geom_def['parameters'][2:5]
+        if any( i >= j for i, j in zip([a, b, c], point)):
+            return False
+
+    if geom_type == "cylinder":
+        r, z = geom_def['parameters'][2:4]
+        if z >= point[2] or any(r >= j for j in point[:2]):
+            return False
+
+    return True
+
 def read_build_file(lines, molecules, topology):
     """
     Parses `lines` of itp format and adds the

polyply/src/meta_molecule.py

@@ -203,6 +203,25 @@ def split_residue(self, split_strings):
         self.relabel_and_redo_res_graph(mapping)
         return mapping
 
+    def copy(self, as_view=False):
+        """
+        Overwrites the networkx copy method. This is needed to
+        correcly set all attributes.
+        """
+        if as_view is True:
+            msg = "Polyply currently does not implement copy as nx graph-views."
+            raise NotImplementedError(msg)
+        G = self.__class__(force_field=self.force_field,
+                           mol_name=self.mol_name)
+        G.molecule = self.molecule
+        G.graph.update(self.graph)
+        G.add_nodes_from((n, d.copy()) for n, d in self._node.items())
+        G.add_edges_from(
+            (u, v, datadict.copy())
+            for u, nbrs in self._adj.items()
+            for v, datadict in nbrs.items())
+        return G
+
     @property
     def search_tree(self):
 

polyply/tests/test_build_file_parser.py

@@ -14,6 +14,7 @@
 """
 Test that build files are properly read.
 """
+import logging
 import textwrap
 import pytest
 import numpy as np
@@ -51,6 +52,41 @@ def test_base_parser_geometry(tokens, _type, expected):
             for result_param, expected_param in zip(result[key][2:], expected[key][2:]):
                 assert result_param == expected_param
 
+@pytest.mark.parametrize('tokens, _type', (
+   # for cylinder radius is not box
+   (["PEO", "63", "154", "in", "4", "8", "8", "6", "7"],
+    "cylinder",),
+   # for cylinder radius is equal
+   (["PEO", "63", "154", "in", "6", "8", "8", "6", "7"],
+    "cylinder",),
+   # for cylinder z is not covered
+   (["PEO", "63", "154", "in", "8", "8", "6", "6", "7"],
+    "cylinder",),
+   # for cylinder z is equal
+   (["PEO", "63", "154", "in", "8", "8", "7", "6", "7"],
+    "cylinder",),
+   # for recangle one of the sides is out
+   (["PEO", "0", "10", "out", "11", "0", "13", "1", "2", "3"],
+    "rectangle",),
+   # for recangle one of the sides is equal
+   (["PEO", "0", "10", "out", "1", "10", "10", "1", "2", "3"],
+    "rectangle",),
+   # for sphere radius is to large/small
+   (["PEO", "0", "10", "in", "0", "12", "13", "5"],
+    "sphere",),
+   # for sphere radius is to equal
+   (["PEO", "0", "10", "in", "5", "12", "13", "5"],
+    "sphere",),
+   ))
+def test_base_parser_geometry_warning(caplog, tokens, _type):
+    with caplog.at_level(logging.WARNING):
+        result = polyply.src.build_file_parser.BuildDirector._base_parser_geometry(tokens, _type)
+        for record in caplog.records:
+            assert record.levelname == "WARNING"
+            break
+        else:
+            assert False
+
 @pytest.fixture
 def test_molecule():
     # dummy vermouth force-field
@@ -115,8 +151,9 @@ def test_system():
                                                     "BB")
   NA = MetaMolecule()
   NA.add_monomer(current=0, resname="NA", connections=[])
+
   molecules = [meta_mol_A, meta_mol_A.copy(),
-               meta_mol_B.copy(), NA, NA.copy(),
+               meta_mol_B, NA, NA.copy(),
                NA.copy(), NA.copy()]
   top = Topology(force_field=force_field)
   top.molecules = molecules
@@ -187,44 +224,96 @@ def test_distance_restraints_error(test_system, line):
    ;
    [ cylinder ]
    ; resname start stop  inside-out  x  y  z  r   z
-   ALA   2    4  in  5  5  5  5  5
+   ALA   1    4  in  5  5  5  4 4
    """,
    [0, 1],
-   [0, 1, 2, 3]),
+   [0, 1, 2]),
    # test that a different molecule is tagged
    ("""
    [ molecule ]
    BB  2  3
    [ rectangle ]
    ; resname start stop  inside-out  x  y  z a b c
-   ALA   3    6  in  5  5  5  5  5 5
+   ALA   3    6  in  5  5  5  4  4 4
    """,
    [2],
-   [2, 3, 4, 5]),
+   [2, 3, 4]),
    # test nothing is tagged based on the molname
+   # combo with mol-idx; molname must be part of
+   # the system though
    ("""
    [ molecule ]
-   CC 1 6
+   BB 0 1
    [ sphere ]
    ; resname start stop  inside-out  x  y  z r
-   ALA   2    4  in  5  5  5  5
+   ALA   2    4  in  5  5  5 4
    """,
    [],
    []),
    ))
 def test_parser(test_system, lines, tagged_mols, tagged_nodes):
-   lines = textwrap.dedent(lines).splitlines()
-   ff = vermouth.forcefield.ForceField(name='test_ff')
-   top = Topology(ff)
-   polyply.src.build_file_parser.read_build_file(lines,
+    lines = textwrap.dedent(lines).splitlines()
+    polyply.src.build_file_parser.read_build_file(lines,
                                                  test_system.molecules,
-                                                 top)
-   for idx, mol in enumerate(test_system.molecules):
-       for node in mol.nodes:
-           if "restraints" in mol.nodes[node]:
-               assert node in tagged_nodes
-               assert idx in tagged_mols
+                                                 test_system)
+    for idx in tagged_mols:
+        mol = test_system.molecules[idx]
+        restr_nodes = nx.get_node_attributes(mol, "restraints")
+        assert list(restr_nodes.keys()) == tagged_nodes
 
+@pytest.mark.parametrize('lines', (
+   # resname not found
+   """
+   [ molecule ]
+   ; name from to
+   AA    0  2
+   ;
+   [ cylinder ]
+   ; resname start stop  inside-out  x  y  z  r   z
+   ALA   2    4  in  5  5  5  4  4
+   PEG   1    2  in  5  5  5  4  4
+   """,
+   # resids don't match resnames
+   """
+   [ molecule ]
+   BB  2  3
+   [ rectangle ]
+   ; resname start stop  inside-out  x  y  z a b c
+   GLU   3    6  in  5  5  5  4  4 4
+   """,
+   # test nothing is tagged based on the molname
+   """
+   [ molecule ]
+   BB 2 3
+   [ sphere ]
+   ; resname start stop  inside-out  x  y  z r
+   GLU   5    6  in  5  5  5  4
+   """,
+   ))
+def test_parser_warnings(caplog, test_system, lines):
+    lines = textwrap.dedent(lines).splitlines()
+    with caplog.at_level(logging.WARNING):
+        polyply.src.build_file_parser.read_build_file(lines,
+                                                      test_system.molecules,
+                                                      test_system)
+        for record in caplog.records:
+            assert record.levelname == "WARNING"
+            break
+        else:
+            assert False
+
+def test_molname_error(test_system):
+    lines = """
+    [ molecule ]
+    CC 0 1
+    [ sphere ]
+    ALA   0    2  in  5  5  5  4
+    """
+    lines = textwrap.dedent(lines).splitlines()
+    with pytest.raises(IOError):
+        polyply.src.build_file_parser.read_build_file(lines,
+                                                      test_system.molecules,
+                                                      test_system)
 
 @pytest.mark.parametrize('lines, expected', (
    # basic test
@@ -258,13 +347,10 @@ def test_parser(test_system, lines, tagged_mols, tagged_nodes):
    )))
 def test_persistence_parsers(test_system, lines, expected):
    lines = textwrap.dedent(lines).splitlines()
-   ff = vermouth.forcefield.ForceField(name='test_ff')
-   top = Topology(ff)
    polyply.src.build_file_parser.read_build_file(lines,
                                                  test_system.molecules,
-                                                 top)
-   for ref, new in zip(expected, top.persistences):
-        print(ref, new)
+                                                 test_system)
+   for ref, new in zip(expected, test_system.persistences):
         for info_ref, info_new in zip(ref[:-1], new[:-1]):
             assert info_ref == info_new
         assert all(ref[-1] == new[-1])

polyply/tests/test_logging.py

@@ -22,7 +22,7 @@
    ("sphere",
     {"resname": "ALA", "start": 9, "stop": 12, "parameters":["in", np.array([10.0, 10.0, 10.0]), 5.0]},
     [],
-    "parsing build file: could not find resid {} with resname ALA in molecule AA.",
+    "parsing build file: could not find residue ALA with resids 9.0,10.0,11.0 in molecule AA.",
     range(9, 12)),
    ))
 def test_tag_nodes_logging(caplog, test_molecule, _type, option, expected, warning, idxs):

polyply/tests/test_meta_molecule.py

@@ -243,3 +243,40 @@ def test_unkown_fromat_error():
         MetaMolecule.from_sequence_file(force_field=ff,
                                         file_path=test_path,
                                         mol_name="test")
+
+def test_copy():
+    """
+    Test if a MetaMol copy has all attributes etc.
+    """
+    file_name = TEST_DATA + "/itp/PEO.itp"
+    edges = [(0,1), (1,2)]
+    nodes = [0, 1, 2]
+    attrs = {0: 'PEO', 1: 'PEO', 2: 'PEO'}
+
+    ff = vermouth.forcefield.ForceField(name='test_ff')
+    name = "PEO"
+    meta_mol = MetaMolecule.from_itp(ff, file_name, name)
+
+    copy_mol = meta_mol.copy()
+    assert list(meta_mol.nodes) == list(copy_mol.nodes)
+    assert meta_mol.force_field == copy_mol.force_field
+    assert meta_mol.molecule == copy_mol.molecule
+    for node in meta_mol.nodes:
+        for key, attr in meta_mol.nodes[node].items():
+            assert attr == copy_mol.nodes[node][key]
+
+def test_copy_as_view():
+    """
+    Test if a MetaMol copy as view raises Error.
+    """
+    file_name = TEST_DATA + "/itp/PEO.itp"
+    edges = [(0,1), (1,2)]
+    nodes = [0, 1, 2]
+    attrs = {0: 'PEO', 1: 'PEO', 2: 'PEO'}
+
+    ff = vermouth.forcefield.ForceField(name='test_ff')
+    name = "PEO"
+    meta_mol = MetaMolecule.from_itp(ff, file_name, name)
+
+    with pytest.raises(NotImplementedError):
+        copy_mol = meta_mol.copy(as_view=True)