EmbedMolecule not respecting double-bond stereochemistry in rings

[jasondbiggs]

The SMILES string refers to trans-cyclodecane, but EmbedMolecule returns a conformer with the cis conformation.

m = Chem.MolFromSmiles("C/1=C\CCCCCCCCCC1")
m = Chem.rdmolops.AddHs(m)
AllChem.EmbedMolecule(m, AllChem.ETKDGv2())
drawit(m)

[JLVarjo]

Hi, I stumbled upon this and ended up making a fix, following how the enforceChirality option works. Briefly, the idea is to look for bond direction tags (ENDUPRIGHT/ENDDOWNRIGHT) around the double bonds, and mark one of the relevant torsions to be either 0 or 180 degrees. Later in the embedding cycles it checks that the generated conformers have the marked torsions correct: if not, the conformer is skipped. Anyway, before I make an actual PR out of this, I'm asking if someone involved to the embedding code would take a look if I'm going to correct direction with this, or missing something critical. The commit is available at JLVarjo@178d7a9

Please note that currently the check breaks some other tests (being too strict for those probably), thus I think it's better to make this under a separate option later.

Cheers, Juuso

[jasondbiggs]

@JLVarjo For my purposes I don't use bond direction for double-bond stereo. I always have the double bond's stereo field set to STEREOZ or STEREOE and have the stereoAtoms set accordingly. To me this makes more sense than the bond direction tags. @greglandrum might be able to say which is the canonical way to check defined double bond stereo in an ROMol.

Also I only have this issue with double bonds inside rings - for non-cyclic double bonds the current embedding routines work fine. I tested the code from your branch in my own embedding code and it worked great, see the image below of cis- and trans-cyclodecene. I encourage you to open a pull request so others with more expertise can opine.

I modified the loop in findStereoBonds to look at the stereo atoms and bond stereo tags, and to only fire for double bonds in sufficiently large rings:

void findStereoBonds(const ROMol &mol,
                     DistGeom::VECT_TORSIONSET &torsionSets,
                     const std::map<int, RDGeom::Point3D> *coordMap) {
  auto rings = mol.getRingInfo();

  if(!rings->isInitialized()) {
    return;
  }

  for(const auto &bond : mol.bonds()) {

    const auto minRing = rings->minBondRingSize(bond->getIdx());

    if(bond->getBondType() == Bond::DOUBLE  && (minRing > 7)) {

      double angle;
      switch(bond->getStereo()) {
        case Bond::STEREOZ:
          angle = 0.0;
          break;
        case Bond::STEREOE:
          angle = M_PI;
          break;
        default:
          continue;
      }

      const auto stereoAtoms = bond->getStereoAtoms();

      if(stereoAtoms.size() != 2) {
        continue;
      }
      torsionSets.emplace_back(
          new DistGeom::TorsionSet(stereoAtoms[0], bond->getBeginAtomIdx(), bond->getEndAtomIdx(), stereoAtoms[1], angle)
      );

    } //end of bond is stereobond condition
  } // end of bond loop
}  // end of _findStereoBonds

[JLVarjo]

I think this is fixed with #5967?