problem with reactions involving radicals

[i-lengyel]

Using a simple bond-breaking reaction smiles should generate two radicals; however, the code adds two extra implicit hydrogens (one for each radicals). Example:
[Cv4:1]-[Cv4:2]>>[Cv3:1].[Cv3:2]

a) Feeding the reaction with 'CCCCC' generates methane ethane, propane and butane rather than methyl, ethyl, ... radicals. It means that there is no atom conservation from left to right, which is certainly incorrect in chemistry.
b) Adding explicit hydrogens to the reactant creates molecules that "looks good" with the explicit hydrogens; however, the molecules are still bad as the code adds an implicit hydrogen to each radicals on the right side. Removing explicit hydrogens certainly does not help, as it will result in a). AssignRadicals function does not do anything as the explicit+implicit hydrogens+C neighbors sum up to the normal valence of C.

Is there a way for a simple fix (like not adding implicit hydrogens in a reaction)?

from rdkit.Chem import AllChem as Chem
from rdkit.Chem import rdmolops as mo
m=Chem.MolFromSmiles('CCCCC')
m=Chem.AddHs(m)
rxn=Chem.ReactionFromSmarts('[Cv4:1]-[Cv4:2]>>[Cv3:1].[Cv3:2]')
ps=rxn.RunReactants((m,)) # produces 8 product pairs, which is 2x that should be
mm=ps[0][0] # pick one product, any
mo.AssignRadicals(mm) # does not do anything, same as without
Chem.SanitizeMol(mm) #have to be used, valence requests does not work without
for a in mm.GetAtoms():
print a.GetSymbol(),a.GetImplicitValence(),a.GetExplicitValence()

results in

C 1 3

H 0 1

H 0 1

H 0 1

[greglandrum]

Hmm, this is a tricky one.
Fixing it would either require introducing some kind of new notation into the SMARTS parser or you doing a bit of extra pre-processing. I need to think about this for a bit.

[i-lengyel]

Greg,

Thanks for looking at it. I was thinking about getting around it by restricting implicit hydrogens with SetNoImplicit, but applying it on the reactants does not affect the the products' atoms. Next I tried is to SetNoImplicit(True) on the products before doing anything with the molecules and that seems work. I guess by default the reaction mapping algorithm when it generates new molecules sets this flag to False. For regular reactions, how biochemists think that could be OK, but for elementary reactions it cannot be right as it does not preserve the number of atoms between the left and right side. Anyway, the algorithm below adjust this:

AddHs to the reactants (make every Hs explicit if they are not already)
RunReactants
for all products:
SetNoImplicit(True)
AssignRadicals
Sanitize (don't know if necessary)
if correct smiles output is needed:
RemoveHs
MolToSmiles will produce the right [CH2]CC-like notations
# end if

end for

This looks like work. I will test a few things with the produced molecules if they act correctly, like pattern matching and such. They should as every Hs are explicit now which should be the case anyway as the whole sequence started with molecules with all explicit Hs. If not, then will do the RemoveHs, MolToSmiles, MolFromSmiles sequence to get to "regular" molecules with both explicit and implicit Hs.

Regards, Istvan

[greglandrum]

Istvan,

As you've noticed, the reaction engine does not really care about stoichiometry. It's designed to be flexible for computational studies, not to model reality.

The procedure you outline above sounds like something that would work for what you're trying to do. It's a clever solution to the problem. I would suggest that it is still a good idea to call sanitize (or at least redo the aromaticity and chirality perception) before generating output. Skipping this step could, under some circumstances, lead to incorrect results.

-greg

[greglandrum]

Closing this, because it looks like the question is answered.