8
$\begingroup$

I am trying to create a molecular graph from SMILES representation of molecules on pytorch using rdkit. The following molecule of caffeine is to be modelled.

Here we see the molecule with bond lengths of each bond in it.

enter image description here

I wish to construct a a pytorch graph for this molecule. For this I would need the node_feat, edge_feat and edge_index data of this molecule. The node_feat is number of nodes x features. edge_feat is number of edge x features . But at the moment my code is including hydrogen bonds as well.

How to neglect/eliminate the hydrogen bonds from this code?

Code:

from rdkit.Chem import AllChem as Chem
from rdkit.Chem import rdMolTransforms
from rdkit.Chem.Draw import rdMolDraw2D, rdDepictor, IPythonConsole
rdDepictor.SetPreferCoordGen(True)
from IPython.display import Image

# smiles = data["smiles"][5:6].values
smiles = ['CN1C=NC2=C1C(=O)N(C(=O)N2C)C']

# 3D for lengths
mol3d = Chem.AddHs(Chem.MolFromSmiles(smiles[0]))
Chem.EmbedMolecule(mol3d, randomSeed=42)
Chem.MMFFOptimizeMolecule(mol3d)

bonds = [(x.GetBeginAtomIdx(), x.GetEndAtomIdx()) for x in mol3d.GetBonds()]

conf = mol3d.GetConformer()

b_lengths = []
b_angles = []

for b in bonds:
    bl = rdMolTransforms.GetBondLength(conf,b[0],b[1]) #GetAngleDeg
    b_lengths.append(str(round(bl,2)))


print(len(b_lengths))
print(b_lengths)
print(len(bonds))
print(bonds)

Output:

25
['1.44', '1.37', '1.32', '1.36', '1.36', '1.42', '1.23', '1.4', '1.41', '1.23', '1.4', '1.45', '1.45', '1.37', '1.37', '1.09', '1.09', '1.09', '1.08', '1.09', '1.09', '1.09', '1.09', '1.09', '1.09']
25
[(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 7), (6, 8), (8, 9), (9, 10), (9, 11), (11, 12), (8, 13), (5, 1), (11, 4), (0, 14), (0, 15), (0, 16), (2, 17), (12, 18), (12, 19), (12, 20), (13, 21), (13, 22), (13, 23)]
$\endgroup$
2
  • $\begingroup$ +1 for the well framed question and details. It’s a delight to answer questions like these when all the required details are provided. $\endgroup$ Commented Apr 15, 2023 at 19:42
  • $\begingroup$ If my answer solved your question, could you please click on the tick mark next to it and mark it as the accepted answer. $\endgroup$ Commented May 18, 2023 at 18:29

1 Answer 1

8
$\begingroup$

You can use the Chem.RemoveHs function.

I attach a snippet below on how you can utilise it in your code.

mol3d_without_Hs = Chem.RemoveHs(mol3d)

bonds = [(x.GetBeginAtomIdx(), x.GetEndAtomIdx()) for x in mol3d_without_Hs.GetBonds()]

conf = mol3d_without_Hs.GetConformer()

b_lengths = []
b_angles = []

for b in bonds:
    bl = rdMolTransforms.GetBondLength(conf,b[0],b[1]) #GetAngleDeg
    b_lengths.append(str(round(bl,2)))

for bnd in mol3d_without_Hs.GetBonds():
    bl = rdMolTransforms.GetBondLength(conf,bnd.GetBeginAtomIdx(),bnd.GetEndAtomIdx())
    # set the bond note on the molecule, which we'll draw:
    mol3d_without_Hs.GetBondWithIdx(bnd.GetIdx()).SetProp('bondNote',f'{bl:.2f}')

Now on a notebook you can visualise the mol3d_without_Hs and you can see the Hs are neglected in the bond lengths:

enter image description here

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .