Friday, November 6, 2020

Mapping atoms to reactants in products made with reaction SMARTS


Reaction SMARTS is an RDKit feature that is very useful for generating reactant products pairs for a given reaction. Unfortunately the algorithm changes the atom order between reactants and products, which creates problems one tries to locate the reaction paths using an interpolation-based algorithm such as nudged elastic band (NEB).

Fortunately, RDKit keeps track of the the change in atom order (thanks to my MS student Julius Seumer for the tip!) and it's easy to reorder the atoms:

def reorder_product(product):
reorder_inverse = [int(atom.GetProp('react_atom_idx')) for atom in product.GetAtoms()]
reorder = len(reorder_inverse)*[0]

for i in range(len(reorder_inverse)):
reorder[reorder_inverse[i]] = i

product = Chem.RenumberAtoms(product, reorder)

return product

If the 3D structures are to be used for interpolation it is important to embed the reactant structure before converting it to product. This keeps the "label-chirality" of groups such as CH2 the same in reactant and products.

Here is a demo notebook

This work is licensed under a Creative Commons Attribution 4.0