RDkit and PySmiles results differ on some SMILES strings

Question

I recently started using both pysmiles and RDkit to parse SMILES strings into molecules. However, I sometimes got different results between the two libraries.

For example, on the molecule described by the string OCCn2c(=N)n(CCOc1ccc(Cl)cc1Cl)c3ccccc23, which is parsed using RDkit into the following molecule:

This molecule has some atom such that RDkit's atom.GetNumImplicitHs() gives a value of $0$, while pysmiles hcount attribute for these atoms is $1$.

This also happens for other molecules, and also with the formal charge value I have seen a similar problem.

Why does this happen? How can I know which library is "right"? Or maybe I didn't understand correctly those features (I'm no chemist, so excuse me for that)? Thanks!

For completeness, the code I used to process the SMILES string with each library is shown below:

• RDkit: molecule = RDkit.Chem.MolFromSmiles(smiles)
• pysmiles: molecule = pysmiles.read_smiles(smiles)

Answer 1

It is a bug in pysmiles.

pysmiles has problems with aromatic nitrogen if the SMILES is in aromatic form (lower case) but not if it is kekulized (upper case and double bonds)

If we add the hydrogens to the molecule, the nitrogens on index 3 and 6 have none, because they have already three bonds.

from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
import pysmiles

s1 = 'OCCn2c(=N)n(CCOc1ccc(Cl)cc1Cl)c3ccccc23'  # aromatic
s2 = 'OCCN2C(=N)N(CCOC1=CC=C(Cl)C=C1Cl)C3=CC=CC=C23'  # kekulized


def show_implicit_h(smiles):
    m = Chem.MolFromSmiles(smiles)
    for atom in m.GetAtoms():
        atom.SetProp('atomLabel', str(atom.GetIdx()))
    m = Chem.AddHs(m)
    return Draw.MolToImage(m, size=(300, 300))


show_implicit_h(s1)

show_implicit_h(s2)

Counting the hydrogen on both smiles shows that pysmiles finds hydrogen in the aromatic, but not in the kekulized SMILES. RDKit finds none in both.

Columns: Index, Symbol, hcount pysmiles, hcount RDKit

def count_implicit_h(s):
    m = Chem.MolFromSmiles(s)

    mol = pysmiles.read_smiles(s)
    hc = mol.nodes(data='hcount')

    for n in hc:
        print(n[0],
              m.GetAtomWithIdx(n[0]).GetSymbol(), n[1],
              m.GetAtomWithIdx(n[0]).GetNumImplicitHs())

aromatic form (count_implicit_h(s1)):

0 O 1 1
1 C 2 2
2 C 2 2
3 N 1 0
4 C 0 0
5 N 1 1
6 N 1 0
7 C 2 2
8 C 2 2
9 O 0 0
10 C 0 0
11 C 1 1
12 C 1 1
13 C 0 0
14 Cl 0 0
15 C 1 1
16 C 0 0
17 Cl 0 0
18 C 0 0
19 C 1 1
20 C 1 1
21 C 1 1
22 C 1 1
23 C 0 0

kekulized form (count_implicit_h(s2)):

0 O 1 1
1 C 2 2
2 C 2 2
3 N 0 0
4 C 0 0
5 N 1 1
6 N 0 0
7 C 2 2
8 C 2 2
9 O 0 0
10 C 0 0
11 C 1 1
12 C 1 1
13 C 0 0
14 Cl 0 0
15 C 1 1
16 C 0 0
17 Cl 0 0
18 C 0 0
19 C 1 1
20 C 1 1
21 C 1 1
22 C 1 1
23 C 0 0

Answer 2

Speculating a comparison of the outputs by both programs may help, I let RDKit process to SMILES to yield a representation including the atom indices. Derived from the cookbook, the MWE for the Jupyter notebook was

from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import rdDepictor
from rdkit.Chem.Draw import rdMolDraw2D

from IPython.display import SVG
IPythonConsole.ipython_useSVG=True  


def mol_with_atom_index(mol):
    for atom in mol.GetAtoms():
        atom.SetAtomMapNum(atom.GetIdx())
    return mol


mol = Chem.MolFromSmiles("OCCn2c(=N)n(CCOc1ccc(Cl)cc1Cl)c3ccccc23")
mol = mol_with_atom_index(mol)
mc = Chem.Mol(mol.ToBinary())


drawer = rdMolDraw2D.MolDraw2DSVG(450, 200) 
drawer.DrawMolecule(mc)
drawer.FinishDrawing()

svg = drawer.GetDrawingText()
display(SVG(svg.replace('svg:','')))

In a second step, I requested pysmiles to list for non-H atoms atom index, atom symbol (first two columns) and number of hydrogen bound (trailing column).

from pysmiles import read_smiles

smiles = "OCCn2c(=N)n(CCOc1ccc(Cl)cc1Cl)c3ccccc23"
molecule = read_smiles(smiles)

pysmiles_list = zip(molecule.nodes(data="element"), molecule.nodes(data="hcount"))
for element in pysmiles_list:
    print(element)

By visual inspection, RDKit's annotated formula and pysmiles' listing agree with each other except for the endocyclic nitrogens (N:3 and N:6):

((0, 'O'), (0, 1))
((1, 'C'), (1, 2))
((2, 'C'), (2, 2))
((3, 'N'), (3, 1))
((4, 'C'), (4, 0))
((5, 'N'), (5, 1))
((6, 'N'), (6, 1))
((7, 'C'), (7, 2))
((8, 'C'), (8, 2))
((9, 'O'), (9, 0))
((10, 'C'), (10, 0))
((11, 'C'), (11, 1))
((12, 'C'), (12, 1))
((13, 'C'), (13, 0))
((14, 'Cl'), (14, 0))
((15, 'C'), (15, 1))
((16, 'C'), (16, 0))
((17, 'Cl'), (17, 0))
((18, 'C'), (18, 0))
((19, 'C'), (19, 1))
((20, 'C'), (20, 1))
((21, 'C'), (21, 1))
((22, 'C'), (22, 1))
((23, 'C'), (23, 0))

A similar listing of atom index, atom type, and number of hydrogens bond to the non-H may be achieved iterating with RDKit over the non-H atoms

for atom in mol.GetAtoms():
    print("{:2} {:2} {}".format(atom.GetIdx(), atom.GetSymbol(),
                                atom.GetTotalNumHs()))

to yield

 0 O  1
 1 C  2
 2 C  2
 3 N  0
 4 C  0
 5 N  1
 6 N  0
 7 C  2
 8 C  2
 9 O  0
10 C  0
11 C  1
12 C  1
13 C  0
14 Cl 0
15 C  1
16 C  0
17 Cl 0
18 C  0
19 C  1
20 C  1
21 C  1
22 C  1
23 C  0

---

As outlined in the answer by @rapelpy, there is a problem with pysmiles which the maintainer(s) might rectify in future. The computations in this answer were performed with pysmiles (version 1.0.1, Jun 25, 2020), and RDKit as packaged for Debian 11 (testing) (version 2020.09.4).