test accounting for explicit hydrogen in the smiles string input

polyply/src/big_smile_mol_processor.py

@@ -104,14 +104,31 @@ def edges_from_bonding_descrpt(self):
             node_graph.nodes[edge[1]]['bonding'] = node_bond_list
             self.meta_molecule.molecule.add_edge(edge[0], edge[1], bonding=bonding)
 
+    def replace_unconsumed_bonding_descrpt(self):
+        """
+        We allow multiple bonding descriptors per atom, which
+        however, are not always consumed. In this case the left
+        over bonding descriptors are replaced by hydrogen atoms.
+        """
+        for node in self.meta_molecule.nodes:
+            graph = self.meta_molecule.nodes[node]['graph']
+            bonding = nx.get_node_attributes(graph, "bonding")
+            for node, bondings in bonding.items():
+                attrs = {attr: graph.nodes[node][attr] for attr in ['resname', 'resid']}
+                attrs['element'] = 'H'
+                for new_id in range(1, len(bondings)+1):
+                    new_node = len(self.meta_molecule.molecule.nodes) + 1
+                    graph.add_edge(node, new_node)
+                    attrs['atomname'] = "H" + str(new_id + len(graph.nodes))
+                    graph.nodes[new_node].update(attrs)
+                    self.meta_molecule.molecule.add_edge(node, new_node)
+                    self.meta_molecule.molecule.nodes[new_node].update(attrs)
+
     def parse(self, big_smile_str):
         res_pattern, residues = big_smile_str.split('.')
         self.meta_molecule = res_pattern_to_meta_mol(res_pattern)
         self.force_field = force_field_from_fragments(residues)
         MapToMolecule(self.force_field).run_molecule(self.meta_molecule)
         self.edges_from_bonding_descrpt()
+        self.replace_unconsumed_bonding_descrpt()
         return self.meta_molecule
-
-# ToDo
-# - replace non consumed bonding descrpt by hydrogen
-# - 

polyply/tests/test_big_smile_mol_proc.py

@@ -52,6 +52,15 @@ def test_generate_edge(bonds_source, bonds_target, edge, btypes):
                         [(0, 1), (0, 2), (2, 3), (3, 4), (2, 5), (2, 6), (4, 7),
                          (4, 8), (4, 9), (9, 10), (10, 11), (9, 12), (9, 13),
                          (11, 14), (11, 15), (11, 16), (16, 17)]),
+                        # uncomsumed bonding IDs; note that this is not the same
+                        # molecule as previous test case. Here one of the OH branches
+                        # and replaces an CH2 group with CH-OH
+                        ("{[#OHter][#PEO]|2[#OHter]}.{#PEO=[>][$1]COC[<],#OHter=[$1][O]}",
+                        [('OHter', 'O H'), ('PEO', 'C O C H H H H'),
+                         ('PEO', 'C O C H H H H'), ('OHter', 'O H')],
+                        [(0, 1), (0, 2), (2, 3), (2, 5), (2, 10), (3, 4),
+                         (4, 6), (4, 7), (4, 17), (8, 9), (8, 11), (8, 14),
+                         (8, 18), (9, 10), (10, 12), (10, 13), (14, 15)]),
                         # simple branched sequence
                         ("{[#Hter][#PE]([#PEO][#Hter])[#PE]([#PEO][#Hter])[#Hter]}.{#Hter=[$]H,#PE=[$]CC[$][$],#PEO=[$]COC[$]}",
                         [('Hter', 'H'), ('PE', 'C C H H H'), ('PEO', 'C O C H H H H'), ('Hter', 'H'),
@@ -75,11 +84,11 @@ def test_generate_edge(bonds_source, bonds_target, edge, btypes):
 ))
 def test_def_big_smile_parser(smile, ref_nodes, ref_edges):
     meta_mol = DefBigSmileParser().parse(smile)
+#    nx.draw_networkx(meta_mol.molecule, with_labels=True, labels=nx.get_node_attributes(meta_mol.molecule, 'element'))
+#    plt.show()
     for node, ref in zip(meta_mol.nodes, ref_nodes):
         assert meta_mol.nodes[node]['resname'] ==  ref[0]
         block_graph = meta_mol.nodes[node]['graph']
         elements = list(nx.get_node_attributes(block_graph, 'element').values())
         assert elements == ref[1].split()
-    #nx.draw_networkx(meta_mol.molecule, with_labels=True, labels=nx.get_node_attributes(meta_mol.molecule, 'element'))
-    #plt.show()
     assert sorted(meta_mol.molecule.edges) == sorted(ref_edges)