Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add export to RDKIT to allow for chemdraw style visualization. #1137

Merged
merged 23 commits into from
Oct 20, 2023
Merged
Show file tree
Hide file tree
Changes from 16 commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
cf759f1
fixed freud bond issue
Jun 9, 2023
457b6ad
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 9, 2023
2f35d1a
Merge branch 'main' into findbonds
daico007 Jun 11, 2023
0096908
Update mbuild/compound.py
daico007 Jun 12, 2023
5302b09
added to_rdkit function with bond order
chrisiacovella Jul 29, 2023
ec1a58f
Merge branch 'mosdef-hub:main' into rdkit_vis
chrisiacovella Jul 30, 2023
332205b
updated bond functions to optionally accept/return bond order
chrisiacovella Jul 30, 2023
fd1ed75
added tests for bond order and converting to rdkit
chrisiacovella Jul 30, 2023
127be7d
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jul 31, 2023
fa7046d
fixed typo in test
chrisiacovella Jul 31, 2023
7d75821
Merge branch 'rdkit_vis' of https://github.com/chrisiacovella/mbuild …
chrisiacovella Jul 31, 2023
acba7c3
fixed typo in test that was missed in last fix
chrisiacovella Jul 31, 2023
0e8cbe7
Merge branch 'main' into rdkit_vis
chrisjonesBSU Oct 18, 2023
3adb518
Check value of bond_order and throw error if invalid, fix typos and s…
chrisjonesBSU Oct 18, 2023
457ae50
move rdkit conv code to conversion
chrisjonesBSU Oct 19, 2023
a4b3921
replace self with compound
chrisjonesBSU Oct 19, 2023
fc023c1
update keys of bond order dict
chrisjonesBSU Oct 19, 2023
510f27d
fix dictionary names and keys
chrisjonesBSU Oct 19, 2023
2ef402c
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Oct 19, 2023
c580b00
add new unit tests, check for type in add_bond
chrisjonesBSU Oct 19, 2023
fa6135f
expand doc strings
chrisjonesBSU Oct 19, 2023
4b14683
add link to rdkit getting started page
chrisjonesBSU Oct 19, 2023
0dd8fa6
add unit test for error in rdkit
chrisjonesBSU Oct 19, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 48 additions & 9 deletions mbuild/compound.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -1204,9 +1211,11 @@ 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(())

Expand Down Expand Up @@ -1246,18 +1255,41 @@ 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
----------
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.
Available options include "default", "single", "double",
"triple", "aromatic" or "unspecified"
"""
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"
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
)

def generate_bonds(self, name_a, name_b, dmin, dmax):
"""Add Bonds between all pairs of types a/b within [dmin, dmax].
Expand Down Expand Up @@ -2654,7 +2686,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"
Expand Down Expand Up @@ -3230,6 +3262,10 @@ def from_parmed(self, structure, coords_only=False, infer_hierarchy=True):
infer_hierarchy=infer_hierarchy,
)

def to_rdkit(self):
chrisjonesBSU marked this conversation as resolved.
Show resolved Hide resolved
"""Create an RDKit RWMol from an mBuild Compound."""
return conversion.to_rdkit(self)

def to_parmed(
self,
box=None,
Expand Down Expand Up @@ -3554,9 +3590,12 @@ 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 "
Expand Down
63 changes: 62 additions & 1 deletion mbuild/conversion.py
Original file line number Diff line number Diff line change
Expand Up @@ -884,10 +884,18 @@
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
Expand Down Expand Up @@ -1749,6 +1757,59 @@
return pybelmol


def to_rdkit(compound):
"""Create an RDKit RWMol from an mBuild Compound."""
rdkit = import_("rdkit")

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable rdkit is not used.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doc strings can probably get added here as well. Should add that we're returning Chem.RWmol, and maybe some examples of useful things to do with the mol and references to rdkit.

Possible examples.

from rdkit.Chem import Draw
img = Draw.MolToImage(temp_mol) # draw a 2D structure
Chem.MolToSmiles(temp_mol) # get smiles string
Chem.RemoveHs(temp_mol) # remove Hydrogens
# list of possible rdkit writers
# link to other working with molecules info (https://www.rdkit.org/docs/GettingStartedInPython.html#working-with-molecules)

Copy link
Contributor

@chrisjonesBSU chrisjonesBSU Oct 19, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good idea. I only added the draw example and the link since mBuild has functionality to remove particles and convert a compound to smiles.

from rdkit import Chem
from rdkit.Chem import AllChem

for particle in compound.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(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,
"default": Chem.BondType.SINGLE,
}

for particle in compound.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 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)
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.

Expand Down
47 changes: 47 additions & 0 deletions mbuild/tests/test_compound.py
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,53 @@ 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
chrisjonesBSU marked this conversation as resolved.
Show resolved Hide resolved

# 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

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
Expand Down
Loading