On Aug 27, 2024, at 14:44, Ingvar Lagerstedt <ing...@nextmovesoftware.com>
wrote:
> To me it would make sense if RDKit removed the aromatic flags for any atom
> that is no longer in a ring when deleting an aromatic atom/bond.
> Alternatively remove the aromatic flag on any non-ring atom when attempting
> to kekulize the structure rather than throwing an exception.
As Noel commented, the toolkit can't make that assumption. Different people
will have different reasons for removing an atom. Some will want to remove
multiple atoms, for example.
Here is one function to remove an atom. It converts the molecule to Kekulé
form, updates the hydrogen counts on the neighboring atoms, and then removes
the specified atom:
from rdkit import Chem
m = Chem.MolFromSmiles('c1ccccc1')
mw = Chem.RWMol(m)
Chem.Kekulize(mw, clearAromaticFlags=True)
atom_idx = 0
atom = mw.GetAtomWithIdx(atom_idx)
for bond in atom.GetBonds():
int_bondtype = int(bond.GetBondType())
assert int_bondtype in (1, 2, 3), "unexpected bond type!"
other_atom = bond.GetOtherAtom(atom)
other_atom.SetNumExplicitHs(
other_atom.GetNumExplicitHs() + int_bondtype)
mw.RemoveAtom(atom_idx)
Chem.SanitizeMol(mw)
print(Chem.MolToSmiles(mw))
Bear in mind that there can be multiple possible Kekulé assignments, while
Chem.Kekulize only picks one. In some cases (not this one, of course) you may
need to apply the above removal method to all distinct assignments (for the
given ring system) in order to get all the valid transformed molecules.
A more correct function should also also take care when the int_bondtype == 1
and the other_atom has a chiral tag, and does not already have a hydrogen,
because you may want to preserve the chiral indicator. Something like this may
work (it's copied and pasted from another function, with a few tweaks to match
the above naming scheme, but I haven't tested it).
num_hs = other_atom.GetTotalNumHs()
if (not num_hs) and other_atom.GetChiralTag():
for bond_i, b in enumerate(other_atom.GetBonds()):
if b.GetIdx() == bond_idx:
break
else:
raise AssertionError("Could not find bond")
want_invert = (bond_i % 2 == 0)
if want_invert:
if other_atom.GetChiralTag() == 2:
other_atom.SetChiralTag(Chem.ChiralType.CHI_TETRAHEDRAL_CW)
else:
other_atom.SetChiralTag(Chem.ChiralType.CHI_TETRAHEDRAL_CCW)
Cheers,
Andrew
da...@dalkescientific.com
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
