How to create a molecular graph from smiles representation

Question

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.

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)]

Answer 1

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: