Assume that I have a molecule as below.
smile = 'C(CN)SS'
mol = Chem.MolFromSmiles(smile)
As can be seen, N has two Hydrogen atoms and S has one Hydrogen atom. I would like to replace those Hydrogens with Carbon.
I can find neighbor atoms by the following code. But, it only returns non-Hydrogen atoms. And, I don't know how to replace Hydrogen with Carbon.
for atom in mol.GetAtoms():
if atom.GetAtomicNum() != 6 and atom.GetTotalNumHs() > 0:
print("Neighbors", [x.GetAtomicNum() for x in atom.GetNeighbors()])
If you have any idea, I really appreciate it to guide me.
The reason you cannot see the hydrogens as Atom objects is because here these hydrogens are modelled as implicit hydrogens. The difference between implicit and explicit hydrogens is summarised nicely in this article. This can also be seen using rdkit:
smiles = 'C(CN)SS'
mol = Chem.MolFromSmiles(smiles)
for atom in mol.GetAtoms():
if atom.GetAtomicNum() != 6 and atom.GetTotalNumHs() > 0:
element = atom.GetSymbol()
implicit = atom.GetNumImplicitHs()
explicit = atom.GetNumExplicitHs()
print(element, 'Implicit:', implicit, 'Explicit:', explicit)
>>> N Implicit: 2 Explicit: 0
>>> S Implicit: 1 Explicit: 0
You can set atoms as explicit easily:
mol = Chem.AddHs(mol)
The easiest way to change substructure is to use the function Chem.ReplaceSubstructs
:
match = Chem.MolFromSmarts('[NH2]')
repl = Chem.MolFromSmarts('N(-C)-C')
new_mol = Chem.ReplaceSubstructs(mol, match, repl)
Okay considering you want to change any hydrogen connected to a non-carbon atom into a carbon you could do something like this:
smiles = 'C(CN)SS'
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
edit = Chem.RWMol(mol)
for atom in edit.GetAtoms():
if atom.GetAtomicNum() != 6:
for nbr in atom.GetNeighbors():
if nbr.GetAtomicNum() == 1:
nbr.SetAtomicNum(6)
Chem.SanitizeMol(edit)
mol = Chem.RemoveHs(edit) # Also converts back to standard Mol