Parsing of hypervalent nitrogen compounds from SMILES/InChI is broken

[mliu49]

Before #1206:

• RMG was able to handle hypervalent N
• RDKit could not really handle hypervalent N (it tries to sanitize and create a non-hypervalent resonance structure)
• OpenBabel could handle hypervalent N
• Parsing of hypervalent SMILES/InChI could be done with OpenBabel

After #1206:

• RMG no longer has atomtypes for hypervalent nitrogen
• Parsing of hypervalent SMILES/InChI breaks with OpenBabel; raises AtomTypeError
• RDKit's treatment of hypervalent N is now similar to RMG's treatment, however RDKit cannot handle all cases.

Examples:

• O=N(=O)C
  • RDKit converts this to one single bonded O and one double bonded O with charge separation
  • OpenBabel converts this to two double bonded O's (raises AtomTypeError)
• O=N(C)=N(C)=O
  • RDKit converts this to N-O single bonds with charge separation and N=N double bond
  • OpenBabel converts this to three double bonds (raises AtomTypeError)
• InChI=1S/C4H10N2O2/c1-3-5(7)6(8)4-2/h3-4H2,1-2H3 (InChI of the prior molecule)
  • RDKit is unable to parse
  • OpenBabel gives the same result

Note: RDKitMol conversion doesn't work properly on current master because of an issue with the order of updateCharge and updateLonePairs. This is a straightforward issue that will be fixed in #1262.

What do we do next?

• Should we re-add atomtypes for hypervalent N?
• Should we add a "sanitation" method similar to RDKit to convert hypervalent N to resonance structures with charge separation?
• Should we prioritize RDKit as a backend? (Currently OpenBabel is tried first for most methods, since it was previously more robust in general.)

[alongd]

I agree that SMILES/InChI parsing is problematic for some cases of hypervalent atoms (N and S). However, I'd like to correct the above claim: to the best of my understanding RMG does have correct atomTypes for all physical hypervalent nitrogen (and sulfur).

The structure O=N(=O)C, for example, isn't a physical molecule representation since the N has a valance of 9 (two double bonds, a single bond, no formal charge). Since in RMG we also use atomTypes to confine the feasible molecule generation space (relevant mainly to resonance structures), we do not have an N5dd atomType (though we do have N5ddc with a formal charge, which does obey the octet rule).

In other words, the three examples above aren't valid physical representations of molecules.

Out of curiosity, could you or @cgrambow elaborate on how these SMILES, which are unphysical representations IMO, were created in the first place?

Other than the examples above, we still experience parsing issues even for valid molecules, e.g., O(N=O)[N]#N, see #735.

[mliu49]

Sorry, I wrote out the SMILES manually. As much as they are un-physical, I would argue that they are still "valid" SMILES which should be corrected to the proper resonance form.

@cgrambow was working with InChI's, which does not specify bond orders. In that case, OpenBabel still interprets these species as I described above. RDKit is capable of interpreting the InChI for O=N(=O)C, but is not able to interpret the second species, which I noted.

[alongd]

We could add such atomTypes to be able to parse "wrong" SMILES, then generate resonance structures and filter out the unrepresentative molecules. However, this quest might be endless. E.g., what about structures like [O-][N+][O+][N-3]#[N+2]? (this is not made up... there are many more examples like this) Where should the line of adding unphysical SMILES cross?

I agree, of course, that InChI aren't specific enough, and cause problems when parsing correct molecule representations. This is an intrinsic InChI problem. I guess that pure InChI isn't designed to represent Lewis structures of hypervalent heteroatoms. Perhaps augmented InChI could solve such cases? (BTW, is augmented InChI something invented by RMG, or is it an IUPAC thing?)

Suggestion: Perhaps we should try both RDKit and OB. If the first fails due to an atomType error, we should try the other. If both fail we could try to assign a correct "similar" atomType to the atom that triggered the error (similar in terms of closest mach of bonds and lone pairs), update nearby atoms accordingly, regenerate InChI, and see if both InChI's match.

[alongd]

Thinking about this more, I agree that adding atomTypes is a more practical solution than the approach I suggested of checking similarities, which probably isn't robust.

My only concern is that we'll needlessly generate structures with wrong atomTypes while applying resonance templates. I don't have a good assessment of how much longer it will take to generate these and then filter them out (keep in mind that in some cases the number of structures could get to hundreds or thousands per species even without the "wrong" atomType additions).

I think we should come up with a way to use these "wrong" atomTypes only for parsing SMILES/InChI (and then resonate to find the correct structure/s), and not when generating new structures. Should we add them as a sub-class of AtomTypes that is only functional for parsing?

A partial list of octet rule violators that were deleted (#1206) is:

atomTypes['Ctc' ] = AtomType('Ctc',  generic=['R','R!H','C','Val4'],  specific=[],
                             single=[0,1,2], allDouble=[0], rDouble=[], oDouble=[], sDouble=[], triple=[1], benzene=[0], lonePairs=[0], charge=[-1,+1])
# examples for Ctc: [CH2-]#[NH+]
atomTypes['N3dc'] = AtomType('N3dc', generic=['R','R!H','N','Val5'],  specific=[],
                             single=[0,1,2], allDouble=[1,2], rDouble=[], oDouble=[], sDouble=[], triple=[0], benzene=[0], lonePairs=[1], charge=[-1,+1])
# *Composite atomType; examples for N3dc: [NH-]=[SH+], [N+]=[S-]=O, C=[N-]=[SH+]
atomTypes['N3tdc'] = AtomType('N3tdc',generic=['R','R!H','N','Val5'],  specific=[],
                             single=[0,1], allDouble=[0,1], rDouble=[], oDouble=[], sDouble=[], triple=[1], benzene=[0], lonePairs=[1], charge=[-1])
# *Composite atomType; examples for N3tdc: [NH-]#[S+]
atomTypes['N5d' ] = AtomType('N5d',   generic=['R','R!H','N','Val5'],  specific=[],
                             single=[0,1,2,3], allDouble=[1], rDouble=[], oDouble=[], sDouble=[], triple=[0], benzene=[0], lonePairs=[0], charge=[0])
# examples for N5d: [NH2]=N
atomTypes['N5dd'] = AtomType('N5dd',  generic=['R','R!H','N','Val5'],  specific=[],
                             single=[0,1], allDouble=[2], rDouble=[], oDouble=[], sDouble=[], triple=[0], benzene=[0], lonePairs=[0], charge=[0])
# examples for N5dd: ON(=O)=O
atomTypes['N5td'] = AtomType('N5td',  generic=['R','R!H','N','Val5'],  specific=[],
                             single=[0], allDouble=[1], rDouble=[], oDouble=[], sDouble=[], triple=[1], benzene=[0], lonePairs=[0], charge=[0])
# examples for N5td: N#N=O
atomTypes['N5tdc'] = AtomType('N5tdc', generic=['R','R!H','N','Val5'],  specific=[],
                             single=[0,1,2,3], allDouble=[0,1,2], rDouble=[], oDouble=[], sDouble=[], triple=[1,2], benzene=[0], lonePairs=[0], charge=[-1])
# *Composite atomType; examples for N5tdc: [NH+]#[NH-]=O
atomTypes['O2sn'] = AtomType('O2sn', generic=['R','R!H','O','Val6'],  specific=[],
                             single=[0,1,2,3], allDouble=[0], rDouble=[], oDouble=[], sDouble=[], triple=[0], benzene=[0], lonePairs=[2], charge=[-1])
# examples for O2sn: C=[NH+][OH2-]
atomTypes['O2dc'] = AtomType('O2dc', generic=['R','R!H','O','Val6'],  specific=[],
                             single=[0,1], allDouble=[1], rDouble=[], oDouble=[], sDouble=[], triple=[0], benzene=[0], lonePairs=[2], charge=[-1])
# examples for O2dc: O=[NH+][O-]=C
atomTypes['O2tc'] = AtomType('O2tc', generic=['R','R!H','O','Val6'],  specific=[],
                             single=[0], allDouble=[0], rDouble=[], oDouble=[], sDouble=[], triple=[1], benzene=[0], lonePairs=[2], charge=[-1])
# examples for O2tc: [O-]#[NH+]

We'll need to add sulfur as well.

[mliu49]

Could we try to update atomtypes, but if it fails, try to generate resonance structures and take the first one that has valid atomtypes?

I don't think it's RMG necessarily needs to use these structure internally, we just need to be able to interpret OB and RDKit output.

[alongd]

Sure! Can we generate a Molecule object without passing through atomTypes?

[github-actions]

This issue is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant issue, otherwise it will automatically be closed in 30 days.