Python program that returns elements from a SMILE string

Question

Given a SMILES string for a molecule, I want to know the elements in the molecule, I then send this to basis_set_exchange to get a "custom" basis set for use in Gaussian.

This seems like something that RDKit can do, or PubChem.

I do know that given a molecular formula, like "H2O" molmass can return the elements. So, it would also suffice to have a program that simply takes a SMILES string and returns a molecular formula, then use molmass to get the elements.

However, this seems like such a trivial exercise, it must already be done. I don't want to make my own script for converting a SMILES to its elements just because I don't want to have to worry about corner cases or reinvent a wheel.

Answer 1

RDKit

This is pretty easy to do in RDKit. If you want the molecular formula, you can just use CalcMolFormula():

from rdkit import Chem
from rdkit.Chem.rdMolDescriptors import CalcMolFormula

# e.g. cysteine
mol = Chem.MolFromSmiles("C([C@@H](C(=O)O)N)S")
formula = CalcMolFormula(mol)

It is also very easy to get all of the unique heteroatoms:

from rdkit import Chem

# e.g. cysteine
mol = Chem.MolFromSmiles("C([C@@H](C(=O)O)N)S")
hetatms = {atom.GetSymbol() for atom in mol.GetAtoms()}

Answer 2

Pybel

Pybel is good for this type of task https://open-babel.readthedocs.io/en/latest/UseTheLibrary/Python_PybelAPI.html

import pybel
global __ATOM_LIST__
 __ATOM_LIST__ = \ 
      ['h',  'he',                                                                                                                                                                                                                                                 
      'li', 'be', 'b',  'c',  'n',  'o',  'f',  'ne',
      'na', 'mg', 'al', 'si', 'p',  's',  'cl', 'ar',
      'k',  'ca', 'sc', 'ti', 'v ', 'cr', 'mn', 'fe', 'co', 'ni', 'cu',
      'zn', 'ga', 'ge', 'as', 'se', 'br', 'kr',
      'rb', 'sr', 'y',  'zr', 'nb', 'mo', 'tc', 'ru', 'rh', 'pd', 'ag',
      'cd', 'in', 'sn', 'sb', 'te', 'i',  'xe',
      'cs', 'ba', 'la', 'ce', 'pr', 'nd', 'pm', 'sm', 'eu', 'gd', 'tb', 'dy',
      'ho', 'er', 'tm', 'yb', 'lu', 'hf', 'ta', 'w',  're', 'os', 'ir', 'pt',
      'au', 'hg', 'tl', 'pb', 'bi', 'po', 'at', 'rn',
      'fr', 'ra', 'ac', 'th', 'pa', 'u',  'np', 'pu']

mol = pybel.readstring("smi", smi)
numbers = [atom.atomicnum for atom in mol]

def from_number_to_element(i):
    return __ATOM_LIST__[i-1].capitalize()

names= [ from_number_to_element(i) for i in numbers]