Reading in wavy bonds from V3K mol blocks

[pechersky]

Hi,

I'm having trouble differentiating wavy single bonds (incident on double bonds) from plain single bonds. I am getting NONE for bond direction in both cases, while I would expect UNKNOWN for the wavy and NONE for the plain. What's the right way to read the mol block in to keep the distinction?

from rdkit import Chem

block = """
  ChemDraw01062412492D

  0  0  0     0  0              0 V3000
M  V30 BEGIN CTAB
M  V30 COUNTS 4 3 0 0 0
M  V30 BEGIN ATOM
M  V30 1 C -0.412500 -0.000000 0.000000 0
M  V30 2 C 0.412500 -0.000000 0.000000 0
M  V30 3 C 0.825000 0.714450 0.000000 0
M  V30 4 C -0.825000 -0.714450 0.000000 0
M  V30 END ATOM
M  V30 BEGIN BOND
M  V30 1 2 1 2 CFG=2
M  V30 2 1 2 3 CFG=2
M  V30 3 1 1 4
M  V30 END BOND
M  V30 END CTAB
M  END

"""

mol = Chem.MolFromMolBlock(block)
bond: Chem.Bond
for bond in mol.GetBonds():
    print(bond.GetBondType(), bond.GetBondDir(), bond.GetStereo())
print()

$ python wavy_check.py
DOUBLE NONE STEREOANY
SINGLE NONE STEREONONE
SINGLE NONE STEREONONE

[rvianello]

You can try and restore the original wedging information from the mol block:

In [3]: mol = Chem.MolFromMolBlock(block)

In [4]: for bond in mol.GetBonds():
   ...:     print(bond.GetBondType(), bond.GetBondDir(), bond.GetStereo())
   ...: 
DOUBLE NONE STEREOANY
SINGLE NONE STEREONONE
SINGLE NONE STEREONONE

In [5]: Chem.ReapplyMolBlockWedging(mol)

In [6]: for bond in mol.GetBonds():
   ...:     print(bond.GetBondType(), bond.GetBondDir(), bond.GetStereo())
   ...: 
DOUBLE EITHERDOUBLE STEREOANY
SINGLE UNKNOWN STEREONONE
SINGLE NONE STEREONONE

or, in case you don't want to modify the mol instance, look for similar information that should be available from the _MolFileBondCfg bond property after the V3000 mol block is parsed (_MolFileBondStereo for V2000).

[pechersky]

Is there documentation anywhere why this data is lost on a regular read, and has to be re-annotated? If one MolFromMolBlock(..., sanitize=False), then the double bond keeps the EITHERDOUBLE. So I'm surprised we go EITHERDOUBLE -> NONE -> EITHERDOUBLE.

[pechersky]

With more detail, I'd like to be able to take this wavy bond molecule, remove the wavy bonds, and then label the double bond as E.

And I'd like to log "there was a wavy bond at X, incident on cross double bond at Y. we cleared the cross and wavy on those, and then were able to assign an E".

[pechersky]

Thanks for the insights, I was also able understand how to do this based on looking at some testChirality tests and reading this discussion

        Chem.ReapplyMolBlockWedging(mol)
        [...]
        cross_bond.SetStereo(Chem.rdchem.BondStereo.STEREONONE)
        cross_bond.SetBondDir(Chem.rdchem.BondDir.NONE)
        wavy_bond.SetBondDir(Chem.rdchem.BondDir.NONE)
        # if we don't clear this, then rdkit's `Chirality::updateDoubleBondNeighbors`
        # will no-op since in addition to the regular API values,
        # it checks this set property
        wavy_bond.ClearProp("_UnknownStereo")
        # without providing the `GetConformer`, the set neighbor directions relies
        # entirely on bond properties -- but that is exactly what we don't have
        # so we must give the 2D coordinates to help along
        Chem.rdmolops.SetDoubleBondNeighborDirections(mol, mol.GetConformer())
        Chem.rdmolops.AssignStereochemistry(mol, False, True)
        stereo_info = next((e for e in Chem.rdmolops.FindPotentialStereo(mol)
                if e.type == Chem.rdchem.StereoType.Bond_Double
                and e.centeredOn == cross_bond.GetIdx()
        ))