diff --git a/polyply/src/build_file_parser.py b/polyply/src/build_file_parser.py index 1e961c00..241ce501 100644 --- a/polyply/src/build_file_parser.py +++ b/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 diff --git a/polyply/src/meta_molecule.py b/polyply/src/meta_molecule.py index 31bbc2e4..da39c8d2 100644 --- a/polyply/src/meta_molecule.py +++ b/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): diff --git a/polyply/tests/test_build_file_parser.py b/polyply/tests/test_build_file_parser.py index fb4447d5..a0630c03 100644 --- a/polyply/tests/test_build_file_parser.py +++ b/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]) diff --git a/polyply/tests/test_logging.py b/polyply/tests/test_logging.py index ed657f96..911de5c1 100644 --- a/polyply/tests/test_logging.py +++ b/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): diff --git a/polyply/tests/test_meta_molecule.py b/polyply/tests/test_meta_molecule.py index 492cb69f..9d532053 100644 --- a/polyply/tests/test_meta_molecule.py +++ b/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)