Fix/380

polyply/src/build_file_parser.py

@@ -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):
@@ -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)
 
@@ -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):
@@ -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.)

polyply/src/check_residue_equivalence.py

@@ -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')

polyply/src/generate_templates.py

@@ -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.
@@ -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]
@@ -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)
@@ -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):
         """

polyply/tests/test_build_file_parser.py

@@ -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

polyply/tests/test_generate_templates.py

@@ -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
@@ -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
@@ -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)

polyply/tests/test_load_library.py

@@ -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
@@ -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'
@@ -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

polyply/tests/test_residue_equivalence.py

@@ -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)