Modify bond graph to include independent particle (#1036)

* initiate bond graph when creating particle

* udpate check in is_independent

* Modify Compound.clone and bond_graph.compose to act on ind node

* Better way to prune bond graph when adding new particle

* exclude port during bond_graph modification

* Modify bond graph behavior regarding removal of nodes and edges
Fix tests of affected functions

* relay kwargs to the gmso conversions method

* trimming code and fix typo

* temporarily disable infer_hierarchy for from_gmso, pending gmso release

* reverse changes relate to unreleased features in gmso

Author: daico007

mbuild/bond_graph.py

@@ -66,6 +66,7 @@ def remove_node(self, node):
         for other_node in self.nodes():
             if node in adj[other_node]:
                 self.remove_edge(node, other_node)
+        del adj[node]
 
     def has_node(self, node):
         """Determine whether the graph contains a node."""
@@ -95,10 +96,6 @@ def remove_edge(self, node1, node2):
         if self.has_node(node1) and self.has_node(node2):
             adj[node1].remove(node2)
             adj[node2].remove(node1)
-            if not adj[node1]:
-                del adj[node1]
-            if not adj[node2]:
-                del adj[node2]
         else:
             raise ValueError(
                 "There is no edge between {} and {}".format(node1, node2)
@@ -151,7 +148,8 @@ def compose(self, graph):
         for node, neighbors in graph._adj.items():
             if self.has_node(node):
                 (adj[node].add(neighbor) for neighbor in neighbors)
-            elif neighbors:
+            else:
+                # Add new node even if it has no bond/neighbor
                 adj[node] = neighbors
 
     def subgraph(self, nodes):

mbuild/compound.py

@@ -175,7 +175,9 @@ def __init__(
         self.labels = OrderedDict()
         self.referrers = set()
 
-        self.bond_graph = None
+        self.bond_graph = BondGraph()
+        self.bond_graph.add_node(self)
+
         self.port_particle = port_particle
 
         self._rigid_id = None
@@ -702,11 +704,13 @@ def add(
             self.children.add(new_child)
             new_child.parent = self
 
-            if new_child.bond_graph is not None:
-                if self.root.bond_graph is None:
-                    self.root.bond_graph = new_child.bond_graph
-                else:
-                    self.root.bond_graph.compose(new_child.bond_graph)
+            if new_child.bond_graph is not None and not isinstance(self, Port):
+                # If anything is added at self level, it is no longer a particle
+                # search for self in self.root.bond_graph and remove self
+                if self.root.bond_graph.has_node(self):
+                    self.root.bond_graph.remove_node(self)
+                # Compose bond_graph of new child
+                self.root.bond_graph.compose(new_child.bond_graph)
 
                 new_child.bond_graph = None
 
@@ -856,7 +860,7 @@ def _remove(self, removed_part):
         if removed_part.rigid_id is not None:
             for ancestor in removed_part.ancestors():
                 ancestor._check_if_contains_rigid_bodies = True
-        if self.root.bond_graph and self.root.bond_graph.has_node(removed_part):
+        if self.root.bond_graph.has_node(removed_part):
             for neighbor in self.root.bond_graph.neighbors(removed_part):
                 self.root.remove_bond((removed_part, neighbor))
             self.root.bond_graph.remove_node(removed_part)
@@ -949,11 +953,8 @@ def direct_bonds(self):
                 "The direct_bonds method can only "
                 "be used on compounds at the bottom of their hierarchy."
             )
-        if not self.root.bond_graph:
-            return iter(())
-        elif self.root.bond_graph.has_node(self):
-            for i in self.root.bond_graph._adj[self]:
-                yield i
+        for i in self.root.bond_graph._adj[self]:
+            yield i
 
     def bonds(self):
         """Return all bonds in the Compound and sub-Compounds.
@@ -1381,7 +1382,7 @@ def is_independent(self):
         if not self.parent:
             # This is the very top level, and hence have to be independent
             return True
-        elif not self.root.bond_graph:
+        elif not self.root.bond_graph.edges():
             # If there is no bond in the top level, then everything is independent
             return True
         else:
@@ -2430,6 +2431,7 @@ def from_gmso(self, topology, coords_only=False, infer_hierarchy=True):
             Set preexisting atoms in compound to coordinates given by Topology.
         infer_hierarchy : bool, optional, default=True
             If True, infer compound hierarchy from Topology residue, to be implemented.
+            Pending new GMSO release.
 
         Returns
         -------
@@ -2439,10 +2441,11 @@ def from_gmso(self, topology, coords_only=False, infer_hierarchy=True):
             topology=topology,
             compound=self,
             coords_only=coords_only,
-            infer_hierarchy=infer_hierarchy,
+            # infer_hierarchy=infer_hierarchy,
+            # TO DO: enable this with new release of GMSO
         )
 
-    def to_gmso(self):
+    def to_gmso(self, **kwargs):
         """Create a GMSO Topology from a mBuild Compound.
 
         Parameters
@@ -2848,6 +2851,9 @@ def _clone(self, clone_of=None, root_container=None):
     def _clone_bonds(self, clone_of=None):
         """Clone the bond of the source compound to clone compound."""
         newone = clone_of[self]
+        newone.bond_graph = BondGraph()
+        for particle in self.particles():
+            newone.bond_graph.add_node(clone_of[particle])
         for c1, c2 in self.bonds():
             try:
                 newone.add_bond((clone_of[c1], clone_of[c2]))

mbuild/conversion.py

@@ -433,7 +433,8 @@ def load_file(
             topology=top,
             compound=compound,
             coords_only=coords_only,
-            infer_hierarchy=infer_hierarchy,
+            # infer_hierarchy=infer_hierarchy,
+            # TODO: enable this with new release of GMSO
         )
 
     # Then pybel reader
@@ -877,6 +878,7 @@ def from_gmso(
         Set preexisting atoms in compound to coordinates given by Topology.
     infer_hierarchy : bool, optional, default=True
         If True, infer compound hierarchy from Topology residue, to be implemented.
+        Pending new GMSO release.
 
     Returns
     -------
@@ -906,9 +908,21 @@ def from_gmso(
 
     # Convert gmso Topology to mbuild Compound
     if not compound:
-        return to_mbuild(topology, **kwargs)
+        return to_mbuild(
+            topology,
+            # infer_hierarchy=infer_hierarchy,
+            # TODO: enable this with new release of GMSO
+            **kwargs,
+        )
     else:
-        compound.add(to_mbuild(topology), **kwargs)
+        compound.add(
+            to_mbuild(
+                topology,
+                # infer_hierarchy=infer_hierarchy),
+                # TODP: enable this with new release of GMSO
+                **kwargs,
+            )
+        )
     return compound
 
 
@@ -1619,13 +1633,15 @@ def _iterate_children(compound, nodes, edges, names_only=False):
     return nodes, edges
 
 
-def to_gmso(compound):
+def to_gmso(compound, box=None, **kwargs):
     """Create a GMSO Topology from a mBuild Compound.
 
     Parameters
     ----------
     compound : mb.Compound
         The mb.Compound to be converted.
+    box : mb.Box, optional, default=None
+        The mb.Box to be converted, if different that compound.box
 
     Returns
     -------
@@ -1634,7 +1650,7 @@ def to_gmso(compound):
     """
     from gmso.external.convert_mbuild import from_mbuild
 
-    return from_mbuild(compound)
+    return from_mbuild(compound, box=None, **kwargs)
 
 
 def to_intermol(compound, molecule_types=None):  # pragma: no cover

mbuild/port.py

@@ -40,8 +40,8 @@ class Port(Compound):
 
     def __init__(self, anchor=None, orientation=None, separation=0):
         super(Port, self).__init__(name="Port", port_particle=True)
+        self.bond_graph = None
         self.anchor = anchor
-
         default_direction = np.array([0, 1, 0])
         if orientation is None:
             orientation = [0, 1, 0]

mbuild/tests/test_compound.py

@@ -1134,14 +1134,20 @@ def test_bond_graph(self, ch3):
         )
 
         ch3_nobonds = mb.clone(ch3)
-        for bond in ch3_nobonds.bonds():
+        bonds_list = list(ch3_nobonds.bonds())
+        for bond in bonds_list:
             ch3_nobonds.remove_bond(bond)
         compound.add(ch3_nobonds)
         assert compound.n_bonds == 3
-        assert not any(
+        assert all(
             compound.bond_graph.has_node(particle)
             for particle in ch3_nobonds.particles()
         )
+        assert not any(
+            compound.bond_graph.has_edge(bond[0], bond[1])
+            for bond in bonds_list
+        )
+        assert compound.bond_graph.has_edge
 
         carbons = list(compound.particles_by_name("C"))
         compound.add_bond((carbons[0], carbons[1]))
@@ -1155,7 +1161,7 @@ def test_bond_graph(self, ch3):
         )
 
         compound.remove_bond((carbons[0], carbons[1]))
-        assert not any(
+        assert all(
             compound.bond_graph.has_node(particle)
             for particle in ch3_nobonds.particles()
         )