From cf759f1f215cefc7a94b3dc7d0166763e6dce1c0 Mon Sep 17 00:00:00 2001 From: Christopher Iacovella Date: Thu, 8 Jun 2023 21:49:56 -0700 Subject: [PATCH 01/19] fixed freud bond issue --- mbuild/compound.py | 19 +++++++++++++++---- mbuild/tests/test_compound.py | 27 +++++++++++++++++++++++---- 2 files changed, 38 insertions(+), 8 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index e114f8f42..24cc1689d 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1306,7 +1306,6 @@ def freud_generate_bonds( name_b, dmin, dmax, - exclude_ii=True, ): """Add Bonds between all pairs of types a/b within [dmin, dmax]. @@ -1320,8 +1319,6 @@ def freud_generate_bonds( The minimum distance (in nm) between Particles for considering a bond dmax : float The maximum distance (in nm) between Particles for considering a bond - exclude_ii : bool, optional, default=True - Whether or not to include neighbors with the same index. Notes ----- @@ -1367,7 +1364,21 @@ def freud_generate_bonds( a_indices.append(i) if part.name == name_b: b_indices.append(i) - + # If we are looking to create bonds between the same species + # then the indices added to a_indices and b_indices will be identical. + # In this case we need to make sure that we don't try to bond a particle + # to itself (i.e., excluded_ii = True). + # If we are looking for bonds between two different species, + # the indices we find for a and b will be distinct, with no overlap. + # However, the way the code is structured, we don't actually pass + # freud the indices, but rather a list of particle positions associated with each set of indices. + # As such, the indices that freud sees will be (0, len(a_indices)) and (0, len(b_indices)), even though + # they represent different actually particles. Thus, toget the right behavior we + # must not exclude particles with the same index, and thus exclude_ii = False. + if name_a == name_b: + exclude_ii=True + else: + exclude_ii=False aq = freud.locality.AABBQuery(freud_box, moved_positions[b_indices]) nlist = aq.query( diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index dc6d46444..936f49c4c 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -721,7 +721,7 @@ def test_freud_generated_bonds_periodicity(self, ch3): ch3_clone.box = mb.Box(lengths=[max(bounding_box.lengths) + 1] * 3) ch3_clone.periodicity = (True, True, True) ch3_clone.freud_generate_bonds( - "H", "H", dmin=0.01, dmax=0.2, exclude_ii=True + "H", "H", dmin=0.01, dmax=0.2 ) assert ch3_clone.n_bonds == 3 + 3 @@ -729,7 +729,7 @@ def test_freud_generated_bonds_periodicity(self, ch3): ch3_clone2.box = mb.Box(lengths=[max(bounding_box.lengths) + 1] * 3) ch3_clone2.periodicity = (True, True, False) ch3_clone2.freud_generate_bonds( - "H", "H", dmin=0.01, dmax=0.2, exclude_ii=True + "H", "H", dmin=0.01, dmax=0.2 ) assert ch3_clone2.n_bonds == 3 + 3 @@ -737,15 +737,34 @@ def test_freud_generated_bonds_periodicity(self, ch3): def test_freud_generate_bonds(self, ch3): bounding_box = ch3.get_boundingbox() ch3.box = mb.Box(lengths=[max(bounding_box.lengths) + 1] * 3) - ch3.freud_generate_bonds("H", "H", dmin=0.01, dmax=0.2, exclude_ii=True) + ch3.freud_generate_bonds("H", "H", dmin=0.01, dmax=0.2) assert ch3.n_bonds == 3 + 3 @pytest.mark.skipif(not has_freud, reason="Freud not installed.") def test_freud_generate_bonds_expected(self, ch3): bounding_box = ch3.get_boundingbox() ch3.box = mb.Box(lengths=[max(bounding_box.lengths) + 1] * 3) - ch3.freud_generate_bonds("H", "H", dmin=0.01, dmax=0.1, exclude_ii=True) + ch3.freud_generate_bonds("H", "H", dmin=0.01, dmax=0.1) assert ch3.n_bonds == 3 + + @pytest.mark.skipif(not has_freud, reason="Freud not installed.") + def test_freud_generate_bonds_mixed(self): + carbon_atom = mb.Compound(name='C', element='C') + + grid_pattern = mb.Grid3DPattern(2, 2, 2) + + grid_pattern.scale([0.25, 0.25, 0.25]) + carbon_list = grid_pattern.apply(carbon_atom) + co_system = mb.Compound(carbon_list) + co_system.box = mb.Box([1,1,1]) + for i, child in enumerate(co_system.children): + if i%2 == 0: + child.name='O' + child.element='O' + + co_system.freud_generate_bonds(name_a= "C", name_b="O", + dmin=0.0, dmax=0.16) + assert co_system.n_bonds == 4 def test_remove_from_box(self, ethane): n_ethanes = 5 From 457b6adba0bfbcbab68d0545e6cdcbb7098d6e22 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 9 Jun 2023 04:58:56 +0000 Subject: [PATCH 02/19] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- mbuild/compound.py | 4 ++-- mbuild/tests/test_compound.py | 25 +++++++++++-------------- 2 files changed, 13 insertions(+), 16 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index 24cc1689d..d68a6783a 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1376,9 +1376,9 @@ def freud_generate_bonds( # they represent different actually particles. Thus, toget the right behavior we # must not exclude particles with the same index, and thus exclude_ii = False. if name_a == name_b: - exclude_ii=True + exclude_ii = True else: - exclude_ii=False + exclude_ii = False aq = freud.locality.AABBQuery(freud_box, moved_positions[b_indices]) nlist = aq.query( diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 936f49c4c..9ab12a226 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -720,17 +720,13 @@ def test_freud_generated_bonds_periodicity(self, ch3): ch3_clone = mb.clone(ch3) ch3_clone.box = mb.Box(lengths=[max(bounding_box.lengths) + 1] * 3) ch3_clone.periodicity = (True, True, True) - ch3_clone.freud_generate_bonds( - "H", "H", dmin=0.01, dmax=0.2 - ) + ch3_clone.freud_generate_bonds("H", "H", dmin=0.01, dmax=0.2) assert ch3_clone.n_bonds == 3 + 3 ch3_clone2 = mb.clone(ch3) ch3_clone2.box = mb.Box(lengths=[max(bounding_box.lengths) + 1] * 3) ch3_clone2.periodicity = (True, True, False) - ch3_clone2.freud_generate_bonds( - "H", "H", dmin=0.01, dmax=0.2 - ) + ch3_clone2.freud_generate_bonds("H", "H", dmin=0.01, dmax=0.2) assert ch3_clone2.n_bonds == 3 + 3 @pytest.mark.skipif(not has_freud, reason="Freud not installed.") @@ -746,24 +742,25 @@ def test_freud_generate_bonds_expected(self, ch3): ch3.box = mb.Box(lengths=[max(bounding_box.lengths) + 1] * 3) ch3.freud_generate_bonds("H", "H", dmin=0.01, dmax=0.1) assert ch3.n_bonds == 3 - + @pytest.mark.skipif(not has_freud, reason="Freud not installed.") def test_freud_generate_bonds_mixed(self): - carbon_atom = mb.Compound(name='C', element='C') + carbon_atom = mb.Compound(name="C", element="C") grid_pattern = mb.Grid3DPattern(2, 2, 2) grid_pattern.scale([0.25, 0.25, 0.25]) carbon_list = grid_pattern.apply(carbon_atom) co_system = mb.Compound(carbon_list) - co_system.box = mb.Box([1,1,1]) + co_system.box = mb.Box([1, 1, 1]) for i, child in enumerate(co_system.children): - if i%2 == 0: - child.name='O' - child.element='O' + if i % 2 == 0: + child.name = "O" + child.element = "O" - co_system.freud_generate_bonds(name_a= "C", name_b="O", - dmin=0.0, dmax=0.16) + co_system.freud_generate_bonds( + name_a="C", name_b="O", dmin=0.0, dmax=0.16 + ) assert co_system.n_bonds == 4 def test_remove_from_box(self, ethane): From 009690824ca77da9bc1ca99958358e7b234fd1e4 Mon Sep 17 00:00:00 2001 From: Co Quach <43968221+daico007@users.noreply.github.com> Date: Mon, 12 Jun 2023 10:33:54 -0500 Subject: [PATCH 03/19] Update mbuild/compound.py --- mbuild/compound.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index d68a6783a..40a3e649e 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1373,7 +1373,7 @@ def freud_generate_bonds( # However, the way the code is structured, we don't actually pass # freud the indices, but rather a list of particle positions associated with each set of indices. # As such, the indices that freud sees will be (0, len(a_indices)) and (0, len(b_indices)), even though - # they represent different actually particles. Thus, toget the right behavior we + # they represent different actually particles. Thus, to get the right behavior we # must not exclude particles with the same index, and thus exclude_ii = False. if name_a == name_b: exclude_ii = True From 5302b09cf91adc994cc528167ec1c8dcef6ffb34 Mon Sep 17 00:00:00 2001 From: chrisiacovella Date: Fri, 28 Jul 2023 23:27:03 -0700 Subject: [PATCH 04/19] added to_rdkit function with bond order --- mbuild/compound.py | 40 ++++++++++++++++++++++++++++++++++++---- mbuild/conversion.py | 10 ++++++++-- 2 files changed, 44 insertions(+), 6 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index 40a3e649e..472044728 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1246,7 +1246,7 @@ def n_bonds(self): ) return sum(1 for _ in self.bonds()) - def add_bond(self, particle_pair): + def add_bond(self, particle_pair, bond_order=None): """Add a bond between two Particles. Parameters @@ -1256,8 +1256,9 @@ def add_bond(self, particle_pair): """ if self.root.bond_graph is None: self.root.bond_graph = BondGraph() - - self.root.bond_graph.add_edge(particle_pair[0], particle_pair[1]) + if bond_order is None: + bond_order = 'default' + self.root.bond_graph.add_edge(particle_pair[0], particle_pair[1], bo=bond_order) def generate_bonds(self, name_a, name_b, dmin, dmax): """Add Bonds between all pairs of types a/b within [dmin, dmax]. @@ -3226,7 +3227,38 @@ def from_parmed(self, structure, coords_only=False, infer_hierarchy=True): coords_only=coords_only, infer_hierarchy=infer_hierarchy, ) - + + def to_rdkit(self): + rdkit = import_("rdkit") + from rdkit import Chem + from rdkit.Chem import AllChem + + temp_mol = Chem.RWMol() + p_dict = {particle: i for i, particle in enumerate(self.particles())} + + bo_dict={ + 1.0: Chem.BondType.SINGLE, + 2.0: Chem.BondType.DOUBLE, + 3.0: Chem.BondType.TRIPLE, + 1.5: Chem.BondType.AROMATIC, + 0.0: Chem.BondType.UNSPECIFIED, + 'default': Chem.BondType.SINGLE} + + + + for particle in self.particles(): + temp_atom = Chem.Atom(particle.element.atomic_number) + temp_atom.SetAtomMapNum(p_dict[particle]) + temp_mol.AddAtom(temp_atom) + + for bond in self.bond_graph.edges(data=True): + bond_indices = (p_dict[bond[0]], p_dict[bond[1]]) + temp_mol.AddBond(*bond_indices) + rdkit_bond = temp_mol.GetBondBetweenAtoms(*bond_indices) + rdkit_bond.SetBondType(bo_dict[bond[2]['bo']]) + + return temp_mol + def to_parmed( self, box=None, diff --git a/mbuild/conversion.py b/mbuild/conversion.py index 50e0f2bb8..e953b8325 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -884,10 +884,16 @@ def from_rdkit(rdkit_mol, compound=None, coords_only=False, smiles_seed=0): part_list.append(part) comp.add(part_list) - + bo_dict={ + Chem.BondType.SINGLE:1.0, + Chem.BondType.DOUBLE:2.0, + Chem.BondType.TRIPLE:3.0, + Chem.BondType.AROMATIC:1.5, + Chem.BondType.UNSPECIFIED:0.0} + for bond in mymol.GetBonds(): comp.add_bond( - [comp[bond.GetBeginAtomIdx()], comp[bond.GetEndAtomIdx()]] + [comp[bond.GetBeginAtomIdx()], comp[bond.GetEndAtomIdx()]], bond_order=bo_dict[bond.GetBondType()] ) return comp From 332205b69f8ec5a6e4e30ac5c7d8e59bf2a689c0 Mon Sep 17 00:00:00 2001 From: chrisiacovella Date: Sun, 30 Jul 2023 12:29:22 -0700 Subject: [PATCH 05/19] updated bond functions to optionally accept/return bond order --- mbuild/compound.py | 49 +++++++++++++++++++++++++++++++++++++--------- 1 file changed, 40 insertions(+), 9 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index 5db6029d3..874bd60d0 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1189,9 +1189,16 @@ def direct_bonds(self): for b1, b2 in self.root.bond_graph.edges(self): yield b2 - def bonds(self): + def bonds(self, return_bond_order=False): """Return all bonds in the Compound and sub-Compounds. - + + Parameters + ---------- + return_bond_order : bool, optional, default=False + Return the bond order of the bond as the 3rd argument in the tuple. + bond order is returned as a dictionary with 'bo' as the key. + If bond order is not set, the value will be set to 'default'. + Yields ------ tuple of mb.Compound @@ -1204,9 +1211,9 @@ def bonds(self): """ if self.root.bond_graph: if self.root == self: - return self.root.bond_graph.edges() + return self.root.bond_graph.edges(data=return_bond_order) else: - return self.root.bond_graph.subgraph(self.particles()).edges() + return self.root.bond_graph.subgraph(self.particles()).edges(data=return_bond_order) else: return iter(()) @@ -1253,6 +1260,8 @@ def add_bond(self, particle_pair, bond_order=None): ---------- particle_pair : indexable object, length=2, dtype=mb.Compound The pair of Particles to add a bond between + bond_order : float, optional, default=None + Bond order of the bond. """ if self.root.bond_graph is None: self.root.bond_graph = BondGraph() @@ -2655,7 +2664,7 @@ def _energy_minimize_openbabel( particle._element = element_from_symbol(particle.name) except ElementError: try: - element_from_name(particle.name) + particle._element = element_from_name(particle.name) except ElementError: raise MBuildError( "No element assigned to {}; element could not be" @@ -3234,6 +3243,23 @@ def to_rdkit(self): from rdkit import Chem from rdkit.Chem import AllChem + + for particle in self.particles(): + if particle.element is None: + try: + particle._element = element_from_symbol(particle.name) + except ElementError: + try: + particle._element = element_from_name(particle.name) + except ElementError: + raise MBuildError( + "No element assigned to {}; element could not be" + "inferred from particle name {}. Cannot perform" + "an energy minimization.".format( + particle, particle.name + ) + ) + temp_mol = Chem.RWMol() p_dict = {particle: i for i, particle in enumerate(self.particles())} @@ -3249,10 +3275,14 @@ def to_rdkit(self): for particle in self.particles(): temp_atom = Chem.Atom(particle.element.atomic_number) - temp_atom.SetAtomMapNum(p_dict[particle]) + + # this next line is necessary to prevent rdkit from adding hydrogens + # this will also set the label to be the lement with particle index + temp_atom.SetProp("atomLabel",f'{temp_atom.GetSymbol()}:{p_dict[particle]}') + temp_mol.AddAtom(temp_atom) - for bond in self.bond_graph.edges(data=True): + for bond in self.bonds(return_bond_order=True): bond_indices = (p_dict[bond[0]], p_dict[bond[1]]) temp_mol.AddBond(*bond_indices) rdkit_bond = temp_mol.GetBondBetweenAtoms(*bond_indices) @@ -3584,9 +3614,10 @@ def _clone_bonds(self, clone_of=None): newone.bond_graph = BondGraph() for particle in self.particles(): newone.bond_graph.add_node(clone_of[particle]) - for c1, c2 in self.bonds(): + for c1, c2, data in self.bonds(return_bond_order=True): try: - newone.add_bond((clone_of[c1], clone_of[c2])) + #bond order is added to the data dictionary as 'bo' + newone.add_bond((clone_of[c1], clone_of[c2]), bond_order=data['bo']) except KeyError: raise MBuildError( "Cloning failed. Compound contains bonds to " From fd1ed75105fa067a9694fbd6fd63eab03a303412 Mon Sep 17 00:00:00 2001 From: chrisiacovella Date: Sun, 30 Jul 2023 12:59:23 -0700 Subject: [PATCH 06/19] added tests for bond order and converting to rdkit --- mbuild/tests/test_compound.py | 49 +++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 18ee707b8..05e7e9715 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -115,6 +115,55 @@ def test_direct_bonds(self, methane): for H in methane.particles_by_name("H"): assert H in bond_particles + def test_bond_order(self, methane): + # test the default behavior + bonds = [bond for bond in Methane().bonds(return_bond_order=True)] + assert bonds[0][2]['bo'] == 'default' + + #test setting bond order when adding a bond + carbon = mb.Compound(pos=[0,0,0], name='C') + carbon_clone = mb.clone(carbon) + C2 = mb.Compound([carbon, carbon_clone]) + C2.add_bond([carbon, carbon_clone], bond_order=2.0) + bonds = [bond for bond in C2.bonds(return_bond_order=True)] + + assert bonds[0][2]['bo'] == 2.0 + + # ensure we propagate bond order information when we clone a compound + C2_clone = mb.clone(C2) + bonds = [bond for bond in C2_clone.bonds(return_bond_order=True)] + + assert bonds[0][2]['bo'] == 2.0 + + @pytest.mark.skipif(not has_rdkit, reason="RDKit is not installed") + def test_convert_to_rdkit(self, methane): + # check basic conversion + rdkit = import_("rdkit") + rdmol = Methane().to_rdkit() + assert isinstance(rdmol, rdkit.Chem.rdchem.RWMol) + + from rdkit import Chem + + rdmol = Methane().to_rdkit() + + atoms = [atom for atom in rdmol.GetAtoms()] + bonds = [bond for bond in rdmol.GetBonds()] + + + bo_dict={ + Chem.BondType.SINGLE:1.0, + Chem.BondType.DOUBLE:2.0, + Chem.BondType.TRIPLE:3.0, + Chem.BondType.AROMATIC:1.5, + Chem.BondType.UNSPECIFIED:0.0} + assert len(atoms) == 5 + assert len(bonds) == 4 + + assert bo_dict[bonds[0].GetBondType()] == 1.0 + assert bo_dict[bonds[1].GetBondType()] == 1.0 + assert bo_dict[bonds[2].GetBondType()] == 1.0 + assert bo_dict[bonds[3].GetBondType()] == 1.0 + def test_n_direct_bonds_parent(self, ethane): with pytest.raises(MBuildError): ethane.n_direct_bonds From 127be7d688f49f76877f2edd24255d79eb0657d6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 31 Jul 2023 00:21:16 +0000 Subject: [PATCH 07/19] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- mbuild/compound.py | 46 ++++++++++++++++++++--------------- mbuild/conversion.py | 18 ++++++++------ mbuild/tests/test_compound.py | 36 +++++++++++++-------------- 3 files changed, 54 insertions(+), 46 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index 874bd60d0..eb0095c1a 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1191,14 +1191,14 @@ def direct_bonds(self): def bonds(self, return_bond_order=False): """Return all bonds in the Compound and sub-Compounds. - + Parameters ---------- return_bond_order : bool, optional, default=False Return the bond order of the bond as the 3rd argument in the tuple. bond order is returned as a dictionary with 'bo' as the key. If bond order is not set, the value will be set to 'default'. - + Yields ------ tuple of mb.Compound @@ -1213,7 +1213,9 @@ def bonds(self, return_bond_order=False): if self.root == self: return self.root.bond_graph.edges(data=return_bond_order) else: - return self.root.bond_graph.subgraph(self.particles()).edges(data=return_bond_order) + return self.root.bond_graph.subgraph(self.particles()).edges( + data=return_bond_order + ) else: return iter(()) @@ -1266,8 +1268,10 @@ def add_bond(self, particle_pair, bond_order=None): if self.root.bond_graph is None: self.root.bond_graph = BondGraph() if bond_order is None: - bond_order = 'default' - self.root.bond_graph.add_edge(particle_pair[0], particle_pair[1], bo=bond_order) + bond_order = "default" + self.root.bond_graph.add_edge( + particle_pair[0], particle_pair[1], bo=bond_order + ) def generate_bonds(self, name_a, name_b, dmin, dmax): """Add Bonds between all pairs of types a/b within [dmin, dmax]. @@ -3237,13 +3241,12 @@ def from_parmed(self, structure, coords_only=False, infer_hierarchy=True): coords_only=coords_only, infer_hierarchy=infer_hierarchy, ) - + def to_rdkit(self): rdkit = import_("rdkit") from rdkit import Chem from rdkit.Chem import AllChem - - + for particle in self.particles(): if particle.element is None: try: @@ -3259,37 +3262,38 @@ def to_rdkit(self): particle, particle.name ) ) - + temp_mol = Chem.RWMol() p_dict = {particle: i for i, particle in enumerate(self.particles())} - bo_dict={ + bo_dict = { 1.0: Chem.BondType.SINGLE, 2.0: Chem.BondType.DOUBLE, 3.0: Chem.BondType.TRIPLE, 1.5: Chem.BondType.AROMATIC, 0.0: Chem.BondType.UNSPECIFIED, - 'default': Chem.BondType.SINGLE} + "default": Chem.BondType.SINGLE, + } - - for particle in self.particles(): temp_atom = Chem.Atom(particle.element.atomic_number) # this next line is necessary to prevent rdkit from adding hydrogens # this will also set the label to be the lement with particle index - temp_atom.SetProp("atomLabel",f'{temp_atom.GetSymbol()}:{p_dict[particle]}') + temp_atom.SetProp( + "atomLabel", f"{temp_atom.GetSymbol()}:{p_dict[particle]}" + ) temp_mol.AddAtom(temp_atom) - + for bond in self.bonds(return_bond_order=True): bond_indices = (p_dict[bond[0]], p_dict[bond[1]]) temp_mol.AddBond(*bond_indices) rdkit_bond = temp_mol.GetBondBetweenAtoms(*bond_indices) - rdkit_bond.SetBondType(bo_dict[bond[2]['bo']]) - + rdkit_bond.SetBondType(bo_dict[bond[2]["bo"]]) + return temp_mol - + def to_parmed( self, box=None, @@ -3616,8 +3620,10 @@ def _clone_bonds(self, clone_of=None): newone.bond_graph.add_node(clone_of[particle]) for c1, c2, data in self.bonds(return_bond_order=True): try: - #bond order is added to the data dictionary as 'bo' - newone.add_bond((clone_of[c1], clone_of[c2]), bond_order=data['bo']) + # bond order is added to the data dictionary as 'bo' + newone.add_bond( + (clone_of[c1], clone_of[c2]), bond_order=data["bo"] + ) except KeyError: raise MBuildError( "Cloning failed. Compound contains bonds to " diff --git a/mbuild/conversion.py b/mbuild/conversion.py index 0e9357084..93b020e29 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -884,16 +884,18 @@ def from_rdkit(rdkit_mol, compound=None, coords_only=False, smiles_seed=0): part_list.append(part) comp.add(part_list) - bo_dict={ - Chem.BondType.SINGLE:1.0, - Chem.BondType.DOUBLE:2.0, - Chem.BondType.TRIPLE:3.0, - Chem.BondType.AROMATIC:1.5, - Chem.BondType.UNSPECIFIED:0.0} - + bo_dict = { + Chem.BondType.SINGLE: 1.0, + Chem.BondType.DOUBLE: 2.0, + Chem.BondType.TRIPLE: 3.0, + Chem.BondType.AROMATIC: 1.5, + Chem.BondType.UNSPECIFIED: 0.0, + } + for bond in mymol.GetBonds(): comp.add_bond( - [comp[bond.GetBeginAtomIdx()], comp[bond.GetEndAtomIdx()]], bond_order=bo_dict[bond.GetBondType()] + [comp[bond.GetBeginAtomIdx()], comp[bond.GetEndAtomIdx()]], + bond_order=bo_dict[bond.GetBondType()], ) return comp diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 05e7e9715..d3989d4c0 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -118,30 +118,30 @@ def test_direct_bonds(self, methane): def test_bond_order(self, methane): # test the default behavior bonds = [bond for bond in Methane().bonds(return_bond_order=True)] - assert bonds[0][2]['bo'] == 'default' - - #test setting bond order when adding a bond - carbon = mb.Compound(pos=[0,0,0], name='C') + assert bonds[0][2]["bo"] == "default" + + # test setting bond order when adding a bond + carbon = mb.Compound(pos=[0, 0, 0], name="C") carbon_clone = mb.clone(carbon) C2 = mb.Compound([carbon, carbon_clone]) C2.add_bond([carbon, carbon_clone], bond_order=2.0) bonds = [bond for bond in C2.bonds(return_bond_order=True)] - assert bonds[0][2]['bo'] == 2.0 - + assert bonds[0][2]["bo"] == 2.0 + # ensure we propagate bond order information when we clone a compound C2_clone = mb.clone(C2) bonds = [bond for bond in C2_clone.bonds(return_bond_order=True)] - assert bonds[0][2]['bo'] == 2.0 - + assert bonds[0][2]["bo"] == 2.0 + @pytest.mark.skipif(not has_rdkit, reason="RDKit is not installed") def test_convert_to_rdkit(self, methane): # check basic conversion rdkit = import_("rdkit") rdmol = Methane().to_rdkit() assert isinstance(rdmol, rdkit.Chem.rdchem.RWMol) - + from rdkit import Chem rdmol = Methane().to_rdkit() @@ -149,21 +149,21 @@ def test_convert_to_rdkit(self, methane): atoms = [atom for atom in rdmol.GetAtoms()] bonds = [bond for bond in rdmol.GetBonds()] - - bo_dict={ - Chem.BondType.SINGLE:1.0, - Chem.BondType.DOUBLE:2.0, - Chem.BondType.TRIPLE:3.0, - Chem.BondType.AROMATIC:1.5, - Chem.BondType.UNSPECIFIED:0.0} + bo_dict = { + Chem.BondType.SINGLE: 1.0, + Chem.BondType.DOUBLE: 2.0, + Chem.BondType.TRIPLE: 3.0, + Chem.BondType.AROMATIC: 1.5, + Chem.BondType.UNSPECIFIED: 0.0, + } assert len(atoms) == 5 assert len(bonds) == 4 - + assert bo_dict[bonds[0].GetBondType()] == 1.0 assert bo_dict[bonds[1].GetBondType()] == 1.0 assert bo_dict[bonds[2].GetBondType()] == 1.0 assert bo_dict[bonds[3].GetBondType()] == 1.0 - + def test_n_direct_bonds_parent(self, ethane): with pytest.raises(MBuildError): ethane.n_direct_bonds From fa7046d9e627ffa788d7ac2df9a5b27ac26b256e Mon Sep 17 00:00:00 2001 From: chrisiacovella Date: Mon, 31 Jul 2023 12:33:58 -0700 Subject: [PATCH 08/19] fixed typo in test --- mbuild/tests/test_compound.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 05e7e9715..12413561a 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -139,13 +139,11 @@ def test_bond_order(self, methane): def test_convert_to_rdkit(self, methane): # check basic conversion rdkit = import_("rdkit") - rdmol = Methane().to_rdkit() + rdmol = methane.to_rdkit() assert isinstance(rdmol, rdkit.Chem.rdchem.RWMol) from rdkit import Chem - rdmol = Methane().to_rdkit() - atoms = [atom for atom in rdmol.GetAtoms()] bonds = [bond for bond in rdmol.GetBonds()] From acba7c3036fb26d10fd312f08a96aa3db4c4aeaf Mon Sep 17 00:00:00 2001 From: chrisiacovella Date: Mon, 31 Jul 2023 12:39:32 -0700 Subject: [PATCH 09/19] fixed typo in test that was missed in last fix --- mbuild/tests/test_compound.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index a2ae85bdd..b84a70510 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -117,7 +117,7 @@ def test_direct_bonds(self, methane): def test_bond_order(self, methane): # test the default behavior - bonds = [bond for bond in Methane().bonds(return_bond_order=True)] + bonds = [bond for bond in methane.bonds(return_bond_order=True)] assert bonds[0][2]["bo"] == "default" # test setting bond order when adding a bond From 3adb5185222c8047493e5c6068fc0487745060ae Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Wed, 18 Oct 2023 17:09:45 -0600 Subject: [PATCH 10/19] Check value of bond_order and throw error if invalid, fix typos and string formatting --- mbuild/compound.py | 30 ++++++++++++++++++++++++------ 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index d9603c4d3..a7b740085 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1264,11 +1264,29 @@ def add_bond(self, particle_pair, bond_order=None): The pair of Particles to add a bond between bond_order : float, optional, default=None Bond order of the bond. + Available options include "default", "single", "double", + "triple", "aromatic" or "unspecified" """ if self.root.bond_graph is None: self.root.bond_graph = BondGraph() if bond_order is None: bond_order = "default" + else: + if bond_order.lower() not in [ + "single", + "double", + "triple", + "aromatic", + "unspecified", + ]: + raise ValueError( + "Invalid bond_order given. Available bond orders are: " + "single", + "double", + "triple", + "aromatic", + "unspecified", + ) self.root.bond_graph.add_edge( particle_pair[0], particle_pair[1], bo=bond_order ) @@ -3245,6 +3263,7 @@ def from_parmed(self, structure, coords_only=False, infer_hierarchy=True): ) def to_rdkit(self): + """Create an RDKit RWMol from an mBuild Compound.""" rdkit = import_("rdkit") from rdkit import Chem from rdkit.Chem import AllChem @@ -3258,11 +3277,10 @@ def to_rdkit(self): particle._element = element_from_name(particle.name) except ElementError: raise MBuildError( - "No element assigned to {}; element could not be" - "inferred from particle name {}. Cannot perform" - "an energy minimization.".format( - particle, particle.name - ) + f"No element assigned to {particle};" + "element could not be" + f"inferred from particle name {particle.name}." + " Cannot perform an energy minimization." ) temp_mol = Chem.RWMol() @@ -3281,7 +3299,7 @@ def to_rdkit(self): temp_atom = Chem.Atom(particle.element.atomic_number) # this next line is necessary to prevent rdkit from adding hydrogens - # this will also set the label to be the lement with particle index + # this will also set the label to be the element with particle index temp_atom.SetProp( "atomLabel", f"{temp_atom.GetSymbol()}:{p_dict[particle]}" ) From 457ae5017464332eb4eec63bf30efd6c5e9afefa Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Wed, 18 Oct 2023 18:05:02 -0600 Subject: [PATCH 11/19] move rdkit conv code to conversion --- mbuild/compound.py | 50 +---------------------------------------- mbuild/conversion.py | 53 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 49 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index a7b740085..13e420056 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -3264,55 +3264,7 @@ def from_parmed(self, structure, coords_only=False, infer_hierarchy=True): def to_rdkit(self): """Create an RDKit RWMol from an mBuild Compound.""" - rdkit = import_("rdkit") - from rdkit import Chem - from rdkit.Chem import AllChem - - for particle in self.particles(): - if particle.element is None: - try: - particle._element = element_from_symbol(particle.name) - except ElementError: - try: - particle._element = element_from_name(particle.name) - except ElementError: - raise MBuildError( - f"No element assigned to {particle};" - "element could not be" - f"inferred from particle name {particle.name}." - " Cannot perform an energy minimization." - ) - - temp_mol = Chem.RWMol() - p_dict = {particle: i for i, particle in enumerate(self.particles())} - - bo_dict = { - 1.0: Chem.BondType.SINGLE, - 2.0: Chem.BondType.DOUBLE, - 3.0: Chem.BondType.TRIPLE, - 1.5: Chem.BondType.AROMATIC, - 0.0: Chem.BondType.UNSPECIFIED, - "default": Chem.BondType.SINGLE, - } - - for particle in self.particles(): - temp_atom = Chem.Atom(particle.element.atomic_number) - - # this next line is necessary to prevent rdkit from adding hydrogens - # this will also set the label to be the element with particle index - temp_atom.SetProp( - "atomLabel", f"{temp_atom.GetSymbol()}:{p_dict[particle]}" - ) - - temp_mol.AddAtom(temp_atom) - - for bond in self.bonds(return_bond_order=True): - bond_indices = (p_dict[bond[0]], p_dict[bond[1]]) - temp_mol.AddBond(*bond_indices) - rdkit_bond = temp_mol.GetBondBetweenAtoms(*bond_indices) - rdkit_bond.SetBondType(bo_dict[bond[2]["bo"]]) - - return temp_mol + return conversion.to_rdkit(self) def to_parmed( self, diff --git a/mbuild/conversion.py b/mbuild/conversion.py index 8d06de806..15631b0ce 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -1757,6 +1757,59 @@ def to_pybel( return pybelmol +def to_rdkit(self): + """Create an RDKit RWMol from an mBuild Compound.""" + rdkit = import_("rdkit") + from rdkit import Chem + from rdkit.Chem import AllChem + + for particle in self.particles(): + if particle.element is None: + try: + particle._element = element_from_symbol(particle.name) + except ElementError: + try: + particle._element = element_from_name(particle.name) + except ElementError: + raise MBuildError( + f"No element assigned to {particle};" + "element could not be" + f"inferred from particle name {particle.name}." + " Cannot perform an energy minimization." + ) + + temp_mol = Chem.RWMol() + p_dict = {particle: i for i, particle in enumerate(self.particles())} + + bo_dict = { + 1.0: Chem.BondType.SINGLE, + 2.0: Chem.BondType.DOUBLE, + 3.0: Chem.BondType.TRIPLE, + 1.5: Chem.BondType.AROMATIC, + 0.0: Chem.BondType.UNSPECIFIED, + "default": Chem.BondType.SINGLE, + } + + for particle in self.particles(): + temp_atom = Chem.Atom(particle.element.atomic_number) + + # this next line is necessary to prevent rdkit from adding hydrogens + # this will also set the label to be the element with particle index + temp_atom.SetProp( + "atomLabel", f"{temp_atom.GetSymbol()}:{p_dict[particle]}" + ) + + temp_mol.AddAtom(temp_atom) + + for bond in self.bonds(return_bond_order=True): + bond_indices = (p_dict[bond[0]], p_dict[bond[1]]) + temp_mol.AddBond(*bond_indices) + rdkit_bond = temp_mol.GetBondBetweenAtoms(*bond_indices) + rdkit_bond.SetBondType(bo_dict[bond[2]["bo"]]) + + return temp_mol + + def to_smiles(compound, backend="pybel"): """Create a SMILES string from an mbuild compound. From a4b392102f1ea1b60d025fe0a98b88de2e6e1258 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Wed, 18 Oct 2023 18:06:09 -0600 Subject: [PATCH 12/19] replace self with compound --- mbuild/conversion.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/mbuild/conversion.py b/mbuild/conversion.py index 15631b0ce..f29805d7c 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -1757,13 +1757,13 @@ def to_pybel( return pybelmol -def to_rdkit(self): +def to_rdkit(compound): """Create an RDKit RWMol from an mBuild Compound.""" rdkit = import_("rdkit") from rdkit import Chem from rdkit.Chem import AllChem - for particle in self.particles(): + for particle in compound.particles(): if particle.element is None: try: particle._element = element_from_symbol(particle.name) @@ -1779,7 +1779,7 @@ def to_rdkit(self): ) temp_mol = Chem.RWMol() - p_dict = {particle: i for i, particle in enumerate(self.particles())} + p_dict = {particle: i for i, particle in enumerate(compound.particles())} bo_dict = { 1.0: Chem.BondType.SINGLE, @@ -1790,7 +1790,7 @@ def to_rdkit(self): "default": Chem.BondType.SINGLE, } - for particle in self.particles(): + for particle in compound.particles(): temp_atom = Chem.Atom(particle.element.atomic_number) # this next line is necessary to prevent rdkit from adding hydrogens @@ -1801,7 +1801,7 @@ def to_rdkit(self): temp_mol.AddAtom(temp_atom) - for bond in self.bonds(return_bond_order=True): + for bond in compound.bonds(return_bond_order=True): bond_indices = (p_dict[bond[0]], p_dict[bond[1]]) temp_mol.AddBond(*bond_indices) rdkit_bond = temp_mol.GetBondBetweenAtoms(*bond_indices) From fc023c1c6d387cbda654967f97860d3cf95a8be2 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Wed, 18 Oct 2023 18:12:44 -0600 Subject: [PATCH 13/19] update keys of bond order dict --- mbuild/conversion.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/mbuild/conversion.py b/mbuild/conversion.py index f29805d7c..13deb15c0 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -1782,11 +1782,11 @@ def to_rdkit(compound): p_dict = {particle: i for i, particle in enumerate(compound.particles())} bo_dict = { - 1.0: Chem.BondType.SINGLE, - 2.0: Chem.BondType.DOUBLE, - 3.0: Chem.BondType.TRIPLE, - 1.5: Chem.BondType.AROMATIC, - 0.0: Chem.BondType.UNSPECIFIED, + "single": Chem.BondType.SINGLE, + "double": Chem.BondType.DOUBLE, + "triple": Chem.BondType.TRIPLE, + "aromatic": Chem.BondType.AROMATIC, + "unspecified": Chem.BondType.UNSPECIFIED, "default": Chem.BondType.SINGLE, } From 510f27d754a7fd2b1b0068911d3a231105e1160c Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 19 Oct 2023 10:34:23 -0600 Subject: [PATCH 14/19] fix dictionary names and keys --- mbuild/compound.py | 5 +++-- mbuild/conversion.py | 18 +++++++++--------- mbuild/tests/test_compound.py | 18 +++++++++--------- 3 files changed, 21 insertions(+), 20 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index 13e420056..960566408 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1273,6 +1273,7 @@ def add_bond(self, particle_pair, bond_order=None): bond_order = "default" else: if bond_order.lower() not in [ + "default", "single", "double", "triple", @@ -1288,7 +1289,7 @@ def add_bond(self, particle_pair, bond_order=None): "unspecified", ) self.root.bond_graph.add_edge( - particle_pair[0], particle_pair[1], bo=bond_order + particle_pair[0], particle_pair[1], bond_order=bond_order ) def generate_bonds(self, name_a, name_b, dmin, dmax): @@ -3594,7 +3595,7 @@ def _clone_bonds(self, clone_of=None): try: # bond order is added to the data dictionary as 'bo' newone.add_bond( - (clone_of[c1], clone_of[c2]), bond_order=data["bo"] + (clone_of[c1], clone_of[c2]), bond_order=data["bond_order"] ) except KeyError: raise MBuildError( diff --git a/mbuild/conversion.py b/mbuild/conversion.py index 13deb15c0..f6a11d542 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -884,18 +884,18 @@ def from_rdkit(rdkit_mol, compound=None, coords_only=False, smiles_seed=0): part_list.append(part) comp.add(part_list) - bo_dict = { - Chem.BondType.SINGLE: 1.0, - Chem.BondType.DOUBLE: 2.0, - Chem.BondType.TRIPLE: 3.0, - Chem.BondType.AROMATIC: 1.5, - Chem.BondType.UNSPECIFIED: 0.0, + bond_order_dict = { + Chem.BondType.SINGLE: "single", + Chem.BondType.DOUBLE: "double", + Chem.BondType.TRIPLE: "triple", + Chem.BondType.AROMATIC: "aromatic", + Chem.BondType.UNSPECIFIED: "unspecified", } for bond in mymol.GetBonds(): comp.add_bond( [comp[bond.GetBeginAtomIdx()], comp[bond.GetEndAtomIdx()]], - bond_order=bo_dict[bond.GetBondType()], + bond_order=bond_order_dict[bond.GetBondType()], ) return comp @@ -1781,7 +1781,7 @@ def to_rdkit(compound): temp_mol = Chem.RWMol() p_dict = {particle: i for i, particle in enumerate(compound.particles())} - bo_dict = { + bond_order_dict = { "single": Chem.BondType.SINGLE, "double": Chem.BondType.DOUBLE, "triple": Chem.BondType.TRIPLE, @@ -1805,7 +1805,7 @@ def to_rdkit(compound): bond_indices = (p_dict[bond[0]], p_dict[bond[1]]) temp_mol.AddBond(*bond_indices) rdkit_bond = temp_mol.GetBondBetweenAtoms(*bond_indices) - rdkit_bond.SetBondType(bo_dict[bond[2]["bo"]]) + rdkit_bond.SetBondType(bond_order_dict[bond[2]["bond_order"]]) return temp_mol diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 832e6daa8..77934e56b 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -118,22 +118,22 @@ def test_direct_bonds(self, methane): def test_bond_order(self, methane): # test the default behavior bonds = [bond for bond in methane.bonds(return_bond_order=True)] - assert bonds[0][2]["bo"] == "default" + assert bonds[0][2]["bond_order"] == "default" # test setting bond order when adding a bond carbon = mb.Compound(pos=[0, 0, 0], name="C") carbon_clone = mb.clone(carbon) C2 = mb.Compound([carbon, carbon_clone]) - C2.add_bond([carbon, carbon_clone], bond_order=2.0) + C2.add_bond([carbon, carbon_clone], bond_order="double") bonds = [bond for bond in C2.bonds(return_bond_order=True)] - assert bonds[0][2]["bo"] == 2.0 + assert bonds[0][2]["bond_order"] == "double" # ensure we propagate bond order information when we clone a compound C2_clone = mb.clone(C2) bonds = [bond for bond in C2_clone.bonds(return_bond_order=True)] - assert bonds[0][2]["bo"] == 2.0 + assert bonds[0][2]["bond_order"] == "double" @pytest.mark.skipif(not has_rdkit, reason="RDKit is not installed") def test_convert_to_rdkit(self, methane): @@ -147,7 +147,7 @@ def test_convert_to_rdkit(self, methane): atoms = [atom for atom in rdmol.GetAtoms()] bonds = [bond for bond in rdmol.GetBonds()] - bo_dict = { + bond_order_dict = { Chem.BondType.SINGLE: 1.0, Chem.BondType.DOUBLE: 2.0, Chem.BondType.TRIPLE: 3.0, @@ -157,10 +157,10 @@ def test_convert_to_rdkit(self, methane): assert len(atoms) == 5 assert len(bonds) == 4 - assert bo_dict[bonds[0].GetBondType()] == 1.0 - assert bo_dict[bonds[1].GetBondType()] == 1.0 - assert bo_dict[bonds[2].GetBondType()] == 1.0 - assert bo_dict[bonds[3].GetBondType()] == 1.0 + assert bond_order_dict[bonds[0].GetBondType()] == 1.0 + assert bond_order_dict[bonds[1].GetBondType()] == 1.0 + assert bond_order_dict[bonds[2].GetBondType()] == 1.0 + assert bond_order_dict[bonds[3].GetBondType()] == 1.0 def test_n_direct_bonds_parent(self, ethane): with pytest.raises(MBuildError): From 2ef402c793772c69c0bf320ce5aeef6880c4329f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 19 Oct 2023 16:34:49 +0000 Subject: [PATCH 15/19] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- mbuild/tests/test_compound.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 77934e56b..b7aee8c91 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -127,13 +127,13 @@ def test_bond_order(self, methane): C2.add_bond([carbon, carbon_clone], bond_order="double") bonds = [bond for bond in C2.bonds(return_bond_order=True)] - assert bonds[0][2]["bond_order"] == "double" + assert bonds[0][2]["bond_order"] == "double" # ensure we propagate bond order information when we clone a compound C2_clone = mb.clone(C2) bonds = [bond for bond in C2_clone.bonds(return_bond_order=True)] - assert bonds[0][2]["bond_order"] == "double" + assert bonds[0][2]["bond_order"] == "double" @pytest.mark.skipif(not has_rdkit, reason="RDKit is not installed") def test_convert_to_rdkit(self, methane): From c580b00c8cb265ae05772a09cf79bf9be33b4494 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 19 Oct 2023 11:41:37 -0600 Subject: [PATCH 16/19] add new unit tests, check for type in add_bond --- mbuild/compound.py | 2 +- mbuild/tests/test_compound.py | 21 +++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index 960566408..695653193 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1272,7 +1272,7 @@ def add_bond(self, particle_pair, bond_order=None): if bond_order is None: bond_order = "default" else: - if bond_order.lower() not in [ + if not isinstance(bond_order, str) or bond_order.lower() not in [ "default", "single", "double", diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index b7aee8c91..02718cf84 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -135,6 +135,27 @@ def test_bond_order(self, methane): assert bonds[0][2]["bond_order"] == "double" + @pytest.mark.parametrize( + "bond_order", + ["single", "double", "triple", "aromatic"], + ) + def test_add_bond_order(self, bond_order): + A_bead = mb.Compound(name="A", pos=[0, 0, 0]) + B_bead = mb.Compound(name="B", pos=[0.5, 0.5, 0]) + comp = mb.Compound([A_bead, B_bead]) + comp.add_bond([A_bead, B_bead], bond_order=bond_order) + for bond in comp.bonds(return_bond_order=True): + assert bond[2]["bond_order"] == bond_order + + def test_add_bad_bond_order(self): + A_bead = mb.Compound(name="A", pos=[0, 0, 0]) + B_bead = mb.Compound(name="B", pos=[0.5, 0.5, 0]) + comp = mb.Compound([A_bead, B_bead]) + with pytest.raises(ValueError): + comp.add_bond([A_bead, B_bead], bond_order=1.0) + with pytest.raises(ValueError): + comp.add_bond([A_bead, B_bead], bond_order="quadruple") + @pytest.mark.skipif(not has_rdkit, reason="RDKit is not installed") def test_convert_to_rdkit(self, methane): # check basic conversion From fa6135f0268ee664045d218ae67fab46e30e6af3 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 19 Oct 2023 12:40:42 -0600 Subject: [PATCH 17/19] expand doc strings --- mbuild/compound.py | 21 ++++++++++++++++++++- mbuild/conversion.py | 12 +++++++++++- 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index 695653193..818a6caa9 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -3264,7 +3264,26 @@ def from_parmed(self, structure, coords_only=False, infer_hierarchy=True): ) def to_rdkit(self): - """Create an RDKit RWMol from an mBuild Compound.""" + """Create an RDKit RWMol from an mBuild Compound. + + Returns + ------- + rdkit.Chem.RWmol + + Notes + ----- + Use this method to utilzie rdkit funcitonality. + + Example + ------- + >>> import mbuild + >>> from rdkit.Chem import Draw + + >>> benzene = mb.load("c1cccc1", smiles=True) + >>> benzene_rdkmol = benzene.to_rdkit() + >>> img = Draw.MolToImage(benzene_rdkmol) + + """ return conversion.to_rdkit(self) def to_parmed( diff --git a/mbuild/conversion.py b/mbuild/conversion.py index f6a11d542..f79302487 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -1758,7 +1758,17 @@ def to_pybel( def to_rdkit(compound): - """Create an RDKit RWMol from an mBuild Compound.""" + """Create an RDKit RWMol from an mBuild Compound. + + Parameters + ---------- + compound : mbuild.Compound; required + The compound to convert to a Chem.RWmol + + Returns + ------- + rdkit.Chem.RWmol + """ rdkit = import_("rdkit") from rdkit import Chem from rdkit.Chem import AllChem From 4b14683dbab05359022b22cf68a6f778a4b01ce9 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 19 Oct 2023 12:45:57 -0600 Subject: [PATCH 18/19] add link to rdkit getting started page --- mbuild/compound.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/mbuild/compound.py b/mbuild/compound.py index 818a6caa9..f407c96e5 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -3283,6 +3283,8 @@ def to_rdkit(self): >>> benzene_rdkmol = benzene.to_rdkit() >>> img = Draw.MolToImage(benzene_rdkmol) + See https://www.rdkit.org/docs/GettingStartedInPython.html + """ return conversion.to_rdkit(self) From 0dd8fa609864ee03eddfb5d9547535a9b5f809ef Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 19 Oct 2023 12:55:47 -0600 Subject: [PATCH 19/19] add unit test for error in rdkit --- mbuild/compound.py | 4 ++++ mbuild/tests/test_compound.py | 8 +++++++- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/mbuild/compound.py b/mbuild/compound.py index f407c96e5..2fe64ba94 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -3273,6 +3273,10 @@ def to_rdkit(self): Notes ----- Use this method to utilzie rdkit funcitonality. + This method only works when the mBuild compound + contains accurate element information. As a result, + this method is not compatible with compounds containing + abstract particles (e.g. coarse-grained systems) Example ------- diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 02718cf84..3c6a00519 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -157,7 +157,7 @@ def test_add_bad_bond_order(self): comp.add_bond([A_bead, B_bead], bond_order="quadruple") @pytest.mark.skipif(not has_rdkit, reason="RDKit is not installed") - def test_convert_to_rdkit(self, methane): + def test_to_rdkit(self, methane): # check basic conversion rdkit = import_("rdkit") rdmol = methane.to_rdkit() @@ -183,6 +183,12 @@ def test_convert_to_rdkit(self, methane): assert bond_order_dict[bonds[2].GetBondType()] == 1.0 assert bond_order_dict[bonds[3].GetBondType()] == 1.0 + @pytest.mark.skipif(not has_rdkit, reason="RDKit is not installed") + def test_to_rdkit_no_elements(self): + comp = mb.Compound(name="_A") + with pytest.raises(MBuildError): + comp.to_rdkit() + def test_n_direct_bonds_parent(self, ethane): with pytest.raises(MBuildError): ethane.n_direct_bonds