Skip to content

Commit

Permalink
Merge pull request #381 from marrink-lab/fix/380
Browse files Browse the repository at this point in the history
Fix/380
  • Loading branch information
csbrasnett authored Jul 26, 2024
2 parents 160e970 + d8d3b66 commit 2e970b4
Show file tree
Hide file tree
Showing 7 changed files with 95 additions and 42 deletions.
22 changes: 17 additions & 5 deletions polyply/src/build_file_parser.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ def __init__(self, molecules, topology):
self.persistence_length = {}
self.templates = {}
self.current_template = None
self.resnames_to_hash = {}

@SectionLineParser.section_parser('molecule')
def _molecule(self, line, lineno=0):
Expand Down Expand Up @@ -152,6 +153,7 @@ def _template_atoms(self, line, lineno=0):
node_name, atype = tokens[0], tokens[1]
position = np.array(tokens[2:], dtype=float)
self.current_template.add_node(node_name,
atomname=node_name,
atype=atype,
position=position)

Expand Down Expand Up @@ -200,14 +202,17 @@ def finalize_section(self, previous_section, ended_section):
# if the volume is not defined yet compute the volume, this still
# can be overwritten by an explicit volume directive later
resname = self.current_template.name
if resname not in self.topology.volumes:
self.topology.volumes[resname] = compute_volume(self.current_template,
coords,
self.topology.nonbond_params,)
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(self.current_template,
node_attr='atomname')
if graph_hash not in self.topology.volumes:
self.topology.volumes[graph_hash] = compute_volume(self.current_template,
coords,
self.topology.nonbond_params,)
# internally a template is defined as vectors from the
# center of geometry
mapped_coords = map_from_CoG(coords)
self.templates[resname] = mapped_coords
self.templates[graph_hash] = mapped_coords
self.resnames_to_hash[resname] = graph_hash
self.current_template = None

def finalize(self, lineno=0):
Expand All @@ -229,6 +234,13 @@ def finalize(self, lineno=0):

super().finalize(lineno=lineno)

# if template graphs and volumes are provided
# make sure that volumes are indexed by the hash
for resname, graph_hash in self.resnames_to_hash.items():
if resname in self.topology.volumes:
self.topology.volumes[graph_hash] = self.topology.volumes[resname]
del self.topology.volumes[resname]

@staticmethod
def _tag_nodes(molecule, keyword, option, molname=""):
resids = np.arange(option['start'], option['stop'], 1.)
Expand Down
6 changes: 1 addition & 5 deletions polyply/src/check_residue_equivalence.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,11 +62,7 @@ def group_residues_by_hash(meta_molecule, template_graphs={}):
dict[`:class:nx.Graph`]
keys are the hash of the graph
"""
unique_graphs = {}
for graph in template_graphs.values():
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(graph, node_attr='atomname')
unique_graphs[graph_hash] = graph

unique_graphs = template_graphs
for node in meta_molecule.nodes:
graph = meta_molecule.nodes[node]["graph"]
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(graph, node_attr='atomname')
Expand Down
16 changes: 9 additions & 7 deletions polyply/src/generate_templates.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@
from .topology import replace_defined_interaction
from .linalg_functions import dih
from .check_residue_equivalence import group_residues_by_hash
from tqdm import tqdm
"""
Processor generating coordinates for all residues of a meta_molecule
matching those in the meta_molecule.molecule attribute.
Expand All @@ -36,6 +35,7 @@ def _extract_template_graphs(meta_molecule, template_graphs={}, skip_filter=Fals
resname = meta_molecule.nodes[node]["resname"]
graph = meta_molecule.nodes[node]["graph"]
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(graph, node_attr='atomname')
meta_molecule.nodes[node]["template"] = graph_hash
if resname in template_graphs:
template_graphs[graph_hash] = graph
del template_graphs[resname]
Expand Down Expand Up @@ -336,8 +336,8 @@ class variable.
self.templates
self.volumes
"""
for resname, template_graph in tqdm(template_graphs.items()):
if resname not in self.templates:
for graph_hash, template_graph in template_graphs.items():
if graph_hash not in self.templates:
block = extract_block(meta_molecule.molecule,
template_graph,
self.topology.defines)
Expand All @@ -363,13 +363,15 @@ class variable.
break
else:
opt_counter += 1

if resname not in self.volumes:
self.volumes[resname] = compute_volume(block,
resname = block.nodes[list(block.nodes)[0]]['resname']
if resname in self.volumes:
self.volumes[graph_hash] = self.volumes[resname]
else:
self.volumes[graph_hash] = compute_volume(block,
coords,
self.topology.nonbond_params)
coords = map_from_CoG(coords)
self.templates[resname] = coords
self.templates[graph_hash] = coords

def run_molecule(self, meta_molecule):
"""
Expand Down
29 changes: 21 additions & 8 deletions polyply/tests/test_build_file_parser.py
Original file line number Diff line number Diff line change
Expand Up @@ -403,12 +403,25 @@ def test_template_volume_parsing(test_system, line, names, edges, positions, out
polyply.src.build_file_parser.read_build_file(lines,
test_system,
test_system.molecules)
for mol in test_system.molecules:
assert len(mol.templates) == len(names)
for idx, name in enumerate(names):
template = mol.templates[name]
for node_pos in positions[idx]:
node = node_pos[0]
assert np.all(np.array(node_pos[1:], dtype=float) == template[node])
# verify results parsing based on hashes
idx = 0
for name, edge_list in zip(names, edges):
template_graph = nx.Graph()
template_graph.add_edges_from(edge_list)
atomnames = {node: node for node in template_graph.nodes}
nx.set_node_attributes(template_graph, atomnames, 'atomname')
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(template_graph,
node_attr='atomname')

assert test_system.volumes == out_vol
mol = test_system.molecules[idx]
assert graph_hash in mol.templates

for node_pos in positions[idx]:
template = mol.templates[graph_hash]
node = node_pos[0]
assert np.all(np.array(node_pos[1:], dtype=float) == template[node])

assert graph_hash in test_system.volumes
assert test_system.volumes[graph_hash] == out_vol[name]

idx += 1
43 changes: 31 additions & 12 deletions polyply/tests/test_generate_templates.py
Original file line number Diff line number Diff line change
Expand Up @@ -153,14 +153,22 @@ def test_extract_block():
len(ff.blocks["GLY"].interactions[inter_type]) == len(new_block.interactions[inter_type])

@staticmethod
def test_run_molecule():
@pytest.mark.parametrize('volumes', (
None,
{"PMMA": 0.55},
))
def test_run_molecule(volumes):
top = polyply.src.topology.Topology.from_gmx_topfile(TEST_DATA / "topology_test" / "system.top", "test")
top.gen_pairs()
if volumes:
top.volumes = volumes
top.convert_nonbond_to_sig_eps()
GenerateTemplates(topology=top, skip_filter=False, max_opt=10).run_molecule(top.molecules[0])
graph = top.molecules[0].nodes[0]['graph']
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(graph, node_attr='atomname')
assert graph_hash in top.volumes
if volumes:
assert top.volumes[graph_hash] == volumes['PMMA']
assert graph_hash in top.molecules[0].templates

@staticmethod
Expand Down Expand Up @@ -309,17 +317,24 @@ def test_compute_volume(lines, coords, volume):
assert np.isclose(new_vol, volume, atol=0.000001)


@pytest.mark.parametrize('resnames, gen_template_graphs, skip_filter', (
@pytest.mark.parametrize('resnames, gen_template_graphs, use_resname, skip_filter', (
# two different residues no template_graphs
(['A', 'B', 'A'], [], False),
(['A', 'B', 'A'], [], False, False),
# two different residues no template_graphs
(['A', 'B', 'A'], [], True),
(['A', 'B', 'A'], [], False, True),
# two different residues one template_graphs
(['A', 'B', 'A'], [1], True),
(['A', 'B', 'A'], [1], False, True),
# two different residues one template_graphs
(['A', 'B', 'A'], [1], False),
(['A', 'B', 'A'], [1], False, False),
# here the template is indexed with the resname
# instead of the hash which needs to be cleared
(['A', 'B', 'A'], [1], True, True),
))
def test_extract_template_graphs(example_meta_molecule, resnames, gen_template_graphs, skip_filter):
def test_extract_template_graphs(example_meta_molecule,
resnames,
gen_template_graphs,
use_resname,
skip_filter):
# set the residue names
for resname, node in zip(resnames, example_meta_molecule.nodes):
example_meta_molecule.nodes[node]['resname'] = resname
Expand All @@ -330,15 +345,19 @@ def test_extract_template_graphs(example_meta_molecule, resnames, gen_template_g
for node in gen_template_graphs:
graph = example_meta_molecule.nodes[node]['graph']
nx.set_node_attributes(graph, True, 'template')
template_graphs[example_meta_molecule.nodes[node]['resname']] = graph
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(graph, node_attr='atomname')
if use_resname:
resname = example_meta_molecule.nodes[node]['resname']
template_graphs[resname] = None
else:
template_graphs[graph_hash] = None

# perfrom the grouping
unique_graphs = _extract_template_graphs(example_meta_molecule, template_graphs, skip_filter)

# check the outcome
assert len(unique_graphs) == 2

for graph in template_graphs.values():
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(graph, node_attr='atomname')
templated = list(nx.get_node_attributes(unique_graphs[graph_hash], 'template').values())
assert all(templated)
# assert that all nodes have the template attribute
for node in example_meta_molecule.nodes:
assert example_meta_molecule.nodes[node].get('template', False)
15 changes: 13 additions & 2 deletions polyply/tests/test_load_library.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
import pytest
from pathlib import Path
from contextlib import contextmanager
import networkx as nx
import vermouth
from polyply import TEST_DATA
from polyply.src.logging import LOGGER
Expand Down Expand Up @@ -72,6 +73,16 @@ def test_read_ff_from_files(caplog):

def test_read_build_options_from_files():

# PMMA template edge_list
edges = [('C1', 'C2'), ('C2', 'C3'), ('C2', 'C4'),
('C4', 'O1'), ('C4', 'O2'), ('O2', 'C5')]
g = nx.Graph()
g.add_edges_from(edges)

atomnames = {node: node for node in g.nodes}
nx.set_node_attributes(g, atomnames, 'atomname')
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(g, node_attr='atomname')

topfile = Path('topology_test/system.top')
bldfile = Path('topology_test/test.bld')
lib_name = '2016H66'
Expand All @@ -83,6 +94,6 @@ def test_read_build_options_from_files():
load_build_files(topology, lib_name, user_files)

# check if build files are parsed
assert topology.volumes == {'PMMA': 1.0}
assert topology.volumes[graph_hash] == 1.0
molecule = topology.molecules[0]
assert molecule.templates
assert graph_hash in molecule.templates
6 changes: 3 additions & 3 deletions polyply/tests/test_residue_equivalence.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,17 +28,17 @@ def test_group_by_hash(example_meta_molecule, resnames, gen_template_graphs):
template_graphs = {}
for node in gen_template_graphs:
graph = example_meta_molecule.nodes[node]['graph']
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(graph, node_attr='atomname')
nx.set_node_attributes(graph, True, 'template')
template_graphs[example_meta_molecule.nodes[node]['resname']] = graph
template_graphs[graph_hash] = graph

# perfrom the grouping
unique_graphs = group_residues_by_hash(example_meta_molecule, template_graphs)

# check the outcome
assert len(unique_graphs) == 2

for graph in template_graphs.values():
graph_hash = nx.algorithms.graph_hashing.weisfeiler_lehman_graph_hash(graph, node_attr='atomname')
for graph_hash in template_graphs:
templated = list(nx.get_node_attributes(unique_graphs[graph_hash], 'template').values())
assert all(templated)

Expand Down

0 comments on commit 2e970b4

Please sign in to comment.