From b97ba7126cfb972bd33f2a1315059f5cbf16a2b8 Mon Sep 17 00:00:00 2001 From: Fabian Gruenewald Date: Mon, 7 Aug 2023 01:00:58 +0200 Subject: [PATCH] filter molecule for templates --- polyply/src/backmap.py | 4 +- polyply/src/gen_coords.py | 5 +- polyply/src/generate_templates.py | 68 +++++++++++++++++------- polyply/tests/test_backmap.py | 6 +-- polyply/tests/test_generate_templates.py | 3 +- 5 files changed, 61 insertions(+), 25 deletions(-) diff --git a/polyply/src/backmap.py b/polyply/src/backmap.py index b643c07b..48851b99 100644 --- a/polyply/src/backmap.py +++ b/polyply/src/backmap.py @@ -166,10 +166,12 @@ def _place_init_coords(self, meta_molecule): ---------- meta_molecule: :class:`polyply.src.MetaMolecule` """ + #print(meta_molecule.templates.keys()) + #print(meta_molecule.templates["PEI2"].keys()) built_nodes = [] for node in meta_molecule.nodes: if meta_molecule.nodes[node]["backmap"]: - resname = meta_molecule.nodes[node]["resname"] + resname = meta_molecule.nodes[node]["template"] cg_coord = meta_molecule.nodes[node]["position"] resid = meta_molecule.nodes[node]["resid"] high_res_atoms = meta_molecule.nodes[node]["graph"].nodes diff --git a/polyply/src/gen_coords.py b/polyply/src/gen_coords.py index 54516cf3..2be1922f 100644 --- a/polyply/src/gen_coords.py +++ b/polyply/src/gen_coords.py @@ -237,8 +237,9 @@ def gen_coords(toppath, grid = np.loadtxt(grid) # do a sanity check - LOGGER.info("checking residue integrity", type="step") - check_residue_equivalence(topology) +# LOGGER.info("checking residue integrity", type="step") +# check_residue_equivalence(topology) + # Build polymer structure LOGGER.info("generating templates", type="step") GenerateTemplates(topology=topology, max_opt=10).run_system(topology) diff --git a/polyply/src/generate_templates.py b/polyply/src/generate_templates.py index 4939353c..84375292 100644 --- a/polyply/src/generate_templates.py +++ b/polyply/src/generate_templates.py @@ -21,6 +21,7 @@ radius_of_gyration) from .topology import replace_defined_interaction from .linalg_functions import dih +from tqdm import tqdm """ Processor generating coordinates for all residues of a meta_molecule matching those in the meta_molecule.molecule attribute. @@ -28,6 +29,9 @@ LOGGER = StyleAdapter(get_logger(__name__)) +def _atoms_match(node1, node2): + return node1["atomname"] == node2["atomname"] + def find_atoms(molecule, attr, value): """ Find all nodes of a `vermouth.molecule.Molecule` that have the @@ -235,7 +239,7 @@ def _relabel_interaction_atoms(interaction, mapping): new_interaction = interaction._replace(atoms=new_atoms) return new_interaction -def extract_block(molecule, resname, defines): +def extract_block(molecule, template_graph, defines): """ Given a `vermouth.molecule` and a `resname` extract the information of a block from the @@ -253,8 +257,8 @@ def extract_block(molecule, resname, defines): ------- :class:vermouth.molecule.Block """ - nodes = find_atoms(molecule, "resname", resname) - resid = molecule.nodes[nodes[0]]["resid"] + #nodes = find_atoms(molecule, "resname", resname) + #resid = molecule.nodes[nodes[0]]["resid"] block = vermouth.molecule.Block() # select all nodes with the same first resid and @@ -263,11 +267,10 @@ def extract_block(molecule, resname, defines): # label in the molecule and in the block for # relabeling the interactions mapping = {} - for node in nodes: + for node in template_graph.nodes: attr_dict = molecule.nodes[node] - if attr_dict["resid"] == resid: - block.add_node(attr_dict["atomname"], **attr_dict) - mapping[node] = attr_dict["atomname"] + block.add_node(attr_dict["atomname"], **attr_dict) + mapping[node] = attr_dict["atomname"] for inter_type in molecule.interactions: for interaction in molecule.interactions[inter_type]: @@ -280,14 +283,36 @@ def extract_block(molecule, resname, defines): "virtual_sites2", "virtual_sites3", "virtual_sites4"]: block.make_edges_from_interaction_type(inter_type) - if not nx.is_connected(block): - msg = ('\n Residue {} with id {} consistes of two disconnected parts. ' - 'Make sure all atoms/particles in a residue are connected by bonds,' - ' constraints or virual-sites.') - raise IOError(msg.format(resname, resid)) + # if not nx.is_connected(block): + # resname = + # msg = ('\n Residue {} with id {} consistes of two disconnected parts. ' + # 'Make sure all atoms/particles in a residue are connected by bonds,' + # ' constraints or virual-sites.') + # raise IOError(msg.format(resname, resid)) return block +def group_by_isomorphism(meta_molecule, template_graphs={}): + """ + Extract all unique fragment graphs from meta_molecule + using the full subgraph isomorphism check. + """ + template_graphs = {} + for node in meta_molecule.nodes: + resname = meta_molecule.nodes[node]["resname"] + graph = meta_molecule.nodes[node]["graph"] + if resname in template_graphs and not nx.is_isomorphic(graph, + template_graphs[resname], + node_match=_atoms_match, + ): + template_name = resname + str(len(template_graphs)) + meta_molecule.nodes[node]["template"] = template_name + template_graphs[template_name] = graph + else: + meta_molecule.nodes[node]["template"] = resname + template_graphs[resname] = graph + return template_graphs + class GenerateTemplates(Processor): """ This processor takes a a class:`polyply.src.MetaMolecule` and @@ -303,7 +328,7 @@ def __init__(self, topology, max_opt, *args, **kwargs): self.volumes = self.topology.volumes self.templates = {} - def gen_templates(self, meta_molecule): + def gen_templates(self, meta_molecule, template_graphs): """ Generate blocks for each unique residue by extracting the block information, placing initial coordinates, and geometry @@ -314,17 +339,19 @@ class variable. Parameters ---------- meta_molecule: :class:`polyply.src.meta_molecule.MetaMolecule` + template_graphs: dict + a dict of graphs corresponding to the templates to be generated Updates --------- self.templates self.volumes """ - resnames = set(nx.get_node_attributes(meta_molecule, "resname").values()) - - for resname in resnames: + #print(len(template_graphs)) + for resname, template_graph in tqdm(template_graphs.items()): if resname not in self.templates: - block = extract_block(meta_molecule.molecule, resname, + block = extract_block(meta_molecule.molecule, + template_graph, self.topology.defines) opt_counter = 0 @@ -362,7 +389,12 @@ def run_molecule(self, meta_molecule): and volume attribute. """ if hasattr(meta_molecule, "templates"): + template_graphs = {res: None for res in meta_molecule.templates} + template_graphs = group_by_isomorphism(meta_molecule, template_graphs) self.templates.update(meta_molecule.templates) - self.gen_templates(meta_molecule) + else: + template_graphs = group_by_isomorphism(meta_molecule) + + self.gen_templates(meta_molecule, template_graphs) meta_molecule.templates = self.templates return meta_molecule diff --git a/polyply/tests/test_backmap.py b/polyply/tests/test_backmap.py index 91aa1ba5..1ecaae83 100644 --- a/polyply/tests/test_backmap.py +++ b/polyply/tests/test_backmap.py @@ -43,13 +43,13 @@ def test_backmapping(): meta_molecule = MetaMolecule(graph) meta_molecule.molecule = molecule - nx.set_node_attributes(meta_molecule, {0: {"resname": "test", + nx.set_node_attributes(meta_molecule, {0: {"resname": "test", "template": "test", "position": np.array([0, 0, 0]), "resid": 1, "backmap": True}, - 1: {"resname": "test", + 1: {"resname": "test", "template": "test", "position": np.array([0, 0, 1.0]), "resid": 2, "backmap": True}, - 2: {"resname": "test", + 2: {"resname": "test", "template": "test", "position": np.array([0, 0, 2.0]), "resid": 3, "backmap": False}}) # test if disordered template works diff --git a/polyply/tests/test_generate_templates.py b/polyply/tests/test_generate_templates.py index ff4bd11f..d82e8201 100644 --- a/polyply/tests/test_generate_templates.py +++ b/polyply/tests/test_generate_templates.py @@ -169,7 +169,8 @@ def test_extract_block(): polyply.src.polyply_parser.read_polyply(lines, ff) block = ff.blocks['test'] molecule = block.to_molecule() - new_block = extract_block(molecule, "GLY", {}) + template_graph = ff.blocks['GLY'].to_molecule() + new_block = extract_block(molecule, template_graph, {}) for node in ff.blocks["GLY"]: atomname = ff.blocks["GLY"].nodes[node]["atomname"] assert ff.blocks["GLY"].nodes[node] == new_block.nodes[atomname]