pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-04 21:54:27 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
c02e53d7aa62dff574da7e6b2dcecfd8f81df023
rdkit/Code/GraphMol/Chirality.cpp
Greg Landrum 02c269f036 Fixes #4115 (#4117)
2021-05-07 10:21:48 -04:00

2618 lines
96 KiB
C++
Raw Blame History

//
// Copyright (C) 2004-2021 Greg Landrum and Rational Discovery LLC
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include <GraphMol/RDKitBase.h>
#include <RDGeneral/Ranking.h>
#include <GraphMol/new_canon.h>
#include <GraphMol/QueryOps.h>
#include <RDGeneral/types.h>
#include <sstream>
#include <set>
#include <algorithm>
#include <RDGeneral/utils.h>
#include <RDGeneral/Invariant.h>
#include <RDGeneral/RDLog.h>
#include <boost/dynamic_bitset.hpp>
#include <Geometry/point.h>
#include "Chirality.h"
// #define VERBOSE_CANON 1
namespace RDKit {
namespace {
bool shouldDetectDoubleBondStereo(const Bond *bond) {
const RingInfo *ri = bond->getOwningMol().getRingInfo();
return (!ri->numBondRings(bond->getIdx()) ||
ri->minBondRingSize(bond->getIdx()) > 7);
}
// ----------------------------------- -----------------------------------
// This algorithm is identical to that used in the CombiCode Mol file
// parser (also developed by RD).
//
//
// SUMMARY:
// Derive a chiral code for an atom that has a wedged (or dashed) bond
// drawn to it.
//
// RETURNS:
// The chiral type
//
// CAVEATS:
// This is careful to ensure that the central atom has 4 neighbors and
// only single bonds to it, but that's about it.
//
// NOTE: this isn't careful at all about checking to make sure that
// things actually *should* be chiral. e.g. if the file has a
// 3-coordinate N with a wedged bond, it will make some erroneous
// assumptions about the chirality.
//
// ----------------------------------- -----------------------------------
Atom::ChiralType atomChiralTypeFromBondDir(const ROMol &mol, const Bond *bond,
const Conformer *conf) {
PRECONDITION(bond, "no bond");
PRECONDITION(conf, "no conformer");
Bond::BondDir bondDir = bond->getBondDir();
PRECONDITION(bondDir == Bond::BEGINWEDGE || bondDir == Bond::BEGINDASH,
"bad bond direction");
// NOTE that according to the CT file spec, wedging assigns chirality
// to the atom at the point of the wedge, (atom 1 in the bond).
const Atom *atom = bond->getBeginAtom();
PRECONDITION(atom, "no atom");
// we can't do anything with atoms that have more than 4 neighbors:
if (atom->getDegree() > 4) {
return Atom::CHI_UNSPECIFIED;
}
const Atom *bondAtom = bond->getEndAtom();
Atom::ChiralType res = Atom::CHI_UNSPECIFIED;
INT_LIST neighborBondIndices;
RDGeom::Point3D centerLoc, tmpPt;
centerLoc = conf->getAtomPos(atom->getIdx());
tmpPt = conf->getAtomPos(bondAtom->getIdx());
centerLoc.z = 0.0;
tmpPt.z = 0.0;
RDGeom::Point3D refVect = centerLoc.directionVector(tmpPt);
//----------------------------------------------------------
//
// start by ensuring that all the bonds to neighboring atoms
// are single bonds and collecting a list of neighbor indices:
//
//----------------------------------------------------------
bool hSeen = false;
neighborBondIndices.push_back(bond->getIdx());
if (bondAtom->getAtomicNum() == 1 && bondAtom->getIsotope() == 0) {
hSeen = true;
}
bool allSingle = true;
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(atom);
while (beg != end) {
const Bond *nbrBond = mol[*beg];
if (nbrBond->getBondType() != Bond::SINGLE) {
allSingle = false;
// break;
}
if (nbrBond != bond) {
if ((nbrBond->getOtherAtom(atom)->getAtomicNum() == 1 &&
nbrBond->getOtherAtom(atom)->getIsotope() == 0)) {
hSeen = true;
}
neighborBondIndices.push_back(nbrBond->getIdx());
}
++beg;
}
size_t nNbrs = neighborBondIndices.size();
//----------------------------------------------------------
//
// Return now if there aren't at least 3 non-H bonds to the atom.
// (we can implicitly add a single H to 3 coordinate atoms, but
// we're horked otherwise).
//
//----------------------------------------------------------
if (nNbrs < 3 || (hSeen && nNbrs < 4)) {
return Atom::CHI_UNSPECIFIED;
}
//----------------------------------------------------------
//
// Continue if there are all single bonds or if we're considering
// 4-coordinate P or S
//
//----------------------------------------------------------
if (allSingle || atom->getAtomicNum() == 15 || atom->getAtomicNum() == 16) {
//------------------------------------------------------------
//
// Here we need to figure out the rotation direction between
// the neighbor bonds and the wedged bond:
//
//------------------------------------------------------------
bool isCCW = true;
INT_LIST::const_iterator bondIter = neighborBondIndices.begin();
++bondIter;
auto bond1 = mol.getBondWithIdx(*bondIter);
int oaid = bond1->getOtherAtom(atom)->getIdx();
tmpPt = conf->getAtomPos(oaid);
tmpPt.z = 0;
auto atomVect0 = centerLoc.directionVector(tmpPt);
auto angle01 = refVect.signedAngleTo(atomVect0);
++bondIter;
auto bond2 = mol.getBondWithIdx(*bondIter);
oaid = bond2->getOtherAtom(atom)->getIdx();
tmpPt = conf->getAtomPos(oaid);
tmpPt.z = 0;
auto atomVect1 = centerLoc.directionVector(tmpPt);
auto angle02 = refVect.signedAngleTo(atomVect1);
// order everything so that looking from the top in a CCW direction we
// have 0, 1, 2
bool swappedIt = false;
if (angle01 > angle02) {
// std::cerr << " swap because " << angle01 << " " << angle02 << " "
// << bond1->getIdx() << "->" << bond2->getIdx() << std::endl;
std::swap(angle01, angle02);
std::swap(bond1, bond2);
std::swap(atomVect0, atomVect1);
swappedIt = true;
}
double firstAngle, secondAngle;
// We proceed differently for 3 and 4 coordinate atoms:
double angle12;
if (nNbrs == 4) {
bool flipIt = false;
// grab the angle to the last neighbor:
++bondIter;
auto bond3 = mol.getBondWithIdx(*bondIter);
oaid = bond3->getOtherAtom(atom)->getIdx();
tmpPt = conf->getAtomPos(oaid);
tmpPt.z = 0;
auto atomVect2 = centerLoc.directionVector(tmpPt);
angle12 = refVect.signedAngleTo(atomVect2);
// find the lowest and second-lowest angle and keep track of
// whether or not we have to do a non-cyclic permutation to
// get there:
if (angle01 < angle02) {
if (angle02 < angle12) {
// order is angle01 -> angle02 -> angle12
firstAngle = angle01;
secondAngle = angle02;
} else if (angle01 < angle12) {
// order is angle01 -> angle12 -> angle02
firstAngle = angle01;
secondAngle = angle12;
flipIt = true;
} else {
// order is angle12 -> angle01 -> angle02
firstAngle = angle12;
secondAngle = angle01;
}
} else if (angle01 < angle12) {
// order is angle02 -> angle01 -> angle12
firstAngle = angle02;
secondAngle = angle01;
flipIt = true;
} else {
if (angle02 < angle12) {
// order is angle02 -> angle12 -> angle01
firstAngle = angle02;
secondAngle = angle12;
} else {
// order is angle12 -> angle02 -> angle01
firstAngle = angle12;
secondAngle = angle02;
flipIt = true;
}
}
if (flipIt) {
isCCW = !isCCW;
}
} else {
// it's three coordinate. Things are a bit different here
// because we have to at least kind of figure out where the
// hydrogen might be.
// before getting started with that, use some of the inchi rules
// for contradictory stereochemistry
// (Table 10 in the InChi v1 technical manual)
angle12 = atomVect0.signedAngleTo(atomVect1);
double angle20 = atomVect1.signedAngleTo(refVect);
// to simplify the code below, pick out the directions of the bonds
// if and only if they start at our atom:
auto dir0 = bondDir;
auto dir1 = (bond1->getBeginAtomIdx() == bond->getBeginAtomIdx())
? bond1->getBondDir()
: Bond::NONE;
auto dir2 = (bond2->getBeginAtomIdx() == bond->getBeginAtomIdx())
? bond2->getBondDir()
: Bond::NONE;
// we know bond 0 has the direction set
// this one is never allowed with different directions:
// 0 2
// \ /
// C
// *
// 1
if (angle01 < (M_PI - 1e-3) && angle12 < (M_PI - 1e-3) &&
angle20 < (M_PI - 1e-3)) {
if ((dir1 != Bond::NONE && dir1 != dir0) ||
(dir2 != Bond::NONE && dir2 != dir0)) {
BOOST_LOG(rdWarningLog)
<< "Warning: conflicting stereochemistry at atom "
<< bond->getBeginAtomIdx() << " ignored."
<< " by rule 1a." << std::endl;
return Atom::CHI_UNSPECIFIED;
}
} else {
// otherwise they cannot all be the same
if (dir1 == dir0 && dir1 == dir2) {
BOOST_LOG(rdWarningLog)
<< "Warning: conflicting stereochemistry at atom "
<< bond->getBeginAtomIdx() << " ignored."
<< " by rule 1b." << std::endl;
return Atom::CHI_UNSPECIFIED;
}
// the remaining cases where stereo is allowed are:
// 0 1 0 2
// * \*/
// 1 - C - 2 C for these two bond1 and bond2 must match
// and
// 1 2 1 0
// * \*/
// 2 - C - 0 C for these two bond0 and bond2 must match
// and
// 2 0 2 1
// * \*/
// 0 - C - 1 C for these two bond0 and bond1 must match
//
if (dir1 != Bond::NONE && dir1 != dir0) {
if ((angle01 >= M_PI) || // last two examples
(angle02 > M_PI && dir2 != Bond::NONE &&
dir1 != dir2) || // top two examples
(angle02 <= M_PI && dir2 != Bond::NONE &&
dir0 != dir2)) { // middle two examples
BOOST_LOG(rdWarningLog)
<< "Warning: conflicting stereochemistry at atom "
<< bond->getBeginAtomIdx() << " ignored."
<< " by rule 2a." << std::endl;
return Atom::CHI_UNSPECIFIED;
}
} else if (dir0 == dir2) {
// all bonds are the same and we already removed the "all
// angles less than 180" case above
BOOST_LOG(rdWarningLog)
<< "Warning: conflicting stereochemistry at atom "
<< bond->getBeginAtomIdx() << " ignored."
<< " by rule 2b." << std::endl;
return Atom::CHI_UNSPECIFIED;
}
if (angle01 < angle02) {
firstAngle = angle01;
secondAngle = angle02;
isCCW = true;
} else {
firstAngle = angle02;
secondAngle = angle01;
isCCW = false;
}
if (secondAngle - firstAngle >= (M_PI - 1e-4)) {
// it's a situation like one of these:
//
// 0 1 0 2
// * \*/
// 1 - C - 2 C
//
// In each of these cases, the implicit H is between atoms 1
// and 2, so we need to flip the rotation direction (go
// around the back).
isCCW = !isCCW;
}
}
}
// reverse the rotation direction if the reference is wedged down:
if (bondDir == Bond::BEGINDASH) {
isCCW = !isCCW;
}
if (swappedIt) {
// we swapped the order of the bonds to simplify the analysis above:
isCCW = !isCCW;
}
// ----------------
//
// We now have the rotation direction using mol-file order.
// We need to convert that into the appropriate label for the
// central atom
//
// ----------------
int nSwaps = atom->getPerturbationOrder(neighborBondIndices);
if (nSwaps % 2) {
isCCW = !isCCW;
}
if (isCCW) {
res = Atom::CHI_TETRAHEDRAL_CCW;
} else {
res = Atom::CHI_TETRAHEDRAL_CW;
}
}
return res;
} // namespace RDKit
Bond::BondDir getOppositeBondDir(Bond::BondDir dir) {
PRECONDITION(dir == Bond::ENDDOWNRIGHT || dir == Bond::ENDUPRIGHT,
"bad bond direction");
switch (dir) {
case Bond::ENDDOWNRIGHT:
return Bond::ENDUPRIGHT;
case Bond::ENDUPRIGHT:
return Bond::ENDDOWNRIGHT;
default:
return Bond::NONE;
}
}
void setBondDirRelativeToAtom(Bond *bond, Atom *atom, Bond::BondDir dir,
bool reverse, boost::dynamic_bitset<> &needsDir) {
PRECONDITION(bond, "bad bond");
PRECONDITION(atom, "bad atom");
PRECONDITION(dir == Bond::ENDUPRIGHT || dir == Bond::ENDDOWNRIGHT, "bad dir");
PRECONDITION(atom == bond->getBeginAtom() || atom == bond->getEndAtom(),
"atom doesn't belong to bond");
RDUNUSED_PARAM(needsDir);
if (bond->getBeginAtom() != atom) {
reverse = !reverse;
}
if (reverse) {
dir = (dir == Bond::ENDUPRIGHT ? Bond::ENDDOWNRIGHT : Bond::ENDUPRIGHT);
}
// to ensure maximum compatibility, even when a bond has unknown stereo (set
// explicitly and recorded in _UnknownStereo property), I will still let a
// direction to be computed. You must check the _UnknownStereo property to
// make sure whether this bond is explicitly set to have no direction info.
// This makes sense because the direction info are all derived from
// coordinates, the _UnknownStereo property is like extra metadata to be
// used with the direction info.
bond->setBondDir(dir);
}
bool isLinearArrangement(const RDGeom::Point3D &v1, const RDGeom::Point3D &v2,
double tol = 0.035) { // tolerance of 2 degrees
return fabs(v2.angleTo(v1) - M_PI) < tol;
}
void updateDoubleBondNeighbors(ROMol &mol, Bond *dblBond, const Conformer *conf,
boost::dynamic_bitset<> &needsDir,
std::vector<unsigned int> &singleBondCounts,
const VECT_INT_VECT &singleBondNbrs) {
// we want to deal only with double bonds:
PRECONDITION(dblBond, "bad bond");
PRECONDITION(dblBond->getBondType() == Bond::DOUBLE, "not a double bond");
if (!needsDir[dblBond->getIdx()]) {
return;
}
needsDir.set(dblBond->getIdx(), 0);
#if 0
std::cerr << "**********************\n";
std::cerr << "**********************\n";
std::cerr << "**********************\n";
std::cerr << "UDBN: " << dblBond->getIdx() << " "
<< dblBond->getBeginAtomIdx() << "=" << dblBond->getEndAtomIdx()
<< "\n";
#endif
ROMol::OEDGE_ITER beg, end;
std::vector<Bond *> followupBonds;
Bond *bond1 = nullptr, *obond1 = nullptr;
bool squiggleBondSeen = false;
bool doubleBondSeen = false;
boost::tie(beg, end) = mol.getAtomBonds(dblBond->getBeginAtom());
while (beg != end) {
Bond *tBond = mol[*beg];
if (tBond == dblBond) {
++beg;
continue;
}
if (tBond->getBondType() == Bond::SINGLE ||
tBond->getBondType() == Bond::AROMATIC) {
// prefer bonds that already have their directionality set
// or that are adjacent to more double bonds:
if (!bond1) {
bond1 = tBond;
} else if (needsDir[tBond->getIdx()]) {
if (singleBondCounts[tBond->getIdx()] >
singleBondCounts[bond1->getIdx()]) {
obond1 = bond1;
bond1 = tBond;
} else {
obond1 = tBond;
}
} else {
obond1 = bond1;
bond1 = tBond;
}
} else if (tBond->getBondType() == Bond::DOUBLE) {
doubleBondSeen = true;
}
int explicit_unknown_stereo;
if ((tBond->getBondType() == Bond::SINGLE ||
tBond->getBondType() == Bond::AROMATIC) &&
(tBond->getBondDir() == Bond::UNKNOWN ||
((tBond->getPropIfPresent<int>(common_properties::_UnknownStereo,
explicit_unknown_stereo) &&
explicit_unknown_stereo)))) {
squiggleBondSeen = true;
break;
}
++beg;
}
// Don't do any direction setting if we've seen a squiggle bond, but do mark
// the double bond as a crossed bond and return
if (!bond1 || squiggleBondSeen || doubleBondSeen) {
if (!doubleBondSeen) {
// FIX: This is the fix for #2649, but it will need to be modified once we
// decide to properly handle allenese
dblBond->setBondDir(Bond::EITHERDOUBLE);
}
return;
}
Bond *bond2 = nullptr, *obond2 = nullptr;
boost::tie(beg, end) = mol.getAtomBonds(dblBond->getEndAtom());
while (beg != end) {
Bond *tBond = mol[*beg];
if (tBond == dblBond) {
++beg;
continue;
}
if (tBond->getBondType() == Bond::SINGLE ||
tBond->getBondType() == Bond::AROMATIC) {
if (!bond2) {
bond2 = tBond;
} else if (needsDir[tBond->getIdx()]) {
if (singleBondCounts[tBond->getIdx()] >
singleBondCounts[bond2->getIdx()]) {
obond2 = bond2;
bond2 = tBond;
} else {
obond2 = tBond;
}
} else {
// we already had a bond2 and we don't need to set the direction
// on the new one, so swap.
obond2 = bond2;
bond2 = tBond;
}
} else if (tBond->getBondType() == Bond::DOUBLE) {
doubleBondSeen = true;
}
int explicit_unknown_stereo;
if (tBond->getBondType() == Bond::SINGLE &&
(tBond->getBondDir() == Bond::UNKNOWN ||
((tBond->getPropIfPresent<int>(common_properties::_UnknownStereo,
explicit_unknown_stereo) &&
explicit_unknown_stereo)))) {
squiggleBondSeen = true;
break;
}
++beg;
}
// Don't do any direction setting if we've seen a squiggle bond, but do mark
// the double bond as a crossed bond and return
if (!bond2 || squiggleBondSeen || doubleBondSeen) {
if (!doubleBondSeen) {
// FIX: This is the fix for #2649, but it will need to be modified once we
// decide to properly handle allenese
dblBond->setBondDir(Bond::EITHERDOUBLE);
}
return;
}
CHECK_INVARIANT(bond1 && bond2, "no bonds found");
bool sameTorsionDir = false;
if (conf) {
RDGeom::Point3D beginP = conf->getAtomPos(dblBond->getBeginAtomIdx());
RDGeom::Point3D endP = conf->getAtomPos(dblBond->getEndAtomIdx());
RDGeom::Point3D bond1P =
conf->getAtomPos(bond1->getOtherAtomIdx(dblBond->getBeginAtomIdx()));
RDGeom::Point3D bond2P =
conf->getAtomPos(bond2->getOtherAtomIdx(dblBond->getEndAtomIdx()));
// check for a linear arrangement of atoms on either end:
bool linear = false;
RDGeom::Point3D p1;
RDGeom::Point3D p2;
p1 = bond1P - beginP;
p2 = endP - beginP;
if (isLinearArrangement(p1, p2)) {
if (!obond1) {
linear = true;
} else {
// one of the bonds was linear; what about the other one?
Bond *tBond = bond1;
bond1 = obond1;
obond1 = tBond;
bond1P = conf->getAtomPos(
bond1->getOtherAtomIdx(dblBond->getBeginAtomIdx()));
p1 = bond1P - beginP;
if (isLinearArrangement(p1, p2)) {
linear = true;
}
}
}
if (!linear) {
p1 = bond2P - endP;
p2 = beginP - endP;
if (isLinearArrangement(p1, p2)) {
if (!obond2) {
linear = true;
} else {
Bond *tBond = bond2;
bond2 = obond2;
obond2 = tBond;
bond2P = conf->getAtomPos(
bond2->getOtherAtomIdx(dblBond->getEndAtomIdx()));
p1 = bond2P - beginP;
if (isLinearArrangement(p1, p2)) {
linear = true;
}
}
}
}
if (linear) {
dblBond->setBondDir(Bond::EITHERDOUBLE);
return;
}
double ang = RDGeom::computeDihedralAngle(bond1P, beginP, endP, bond2P);
if (ang < M_PI / 2) {
sameTorsionDir = false;
} else {
sameTorsionDir = true;
}
// std::cerr << " angle: " << ang << " sameTorsionDir: " << sameTorsionDir
// << "\n";
} else {
if (dblBond->getStereo() == Bond::STEREOCIS ||
dblBond->getStereo() == Bond::STEREOZ) {
sameTorsionDir = false;
} else if (dblBond->getStereo() == Bond::STEREOTRANS ||
dblBond->getStereo() == Bond::STEREOE) {
sameTorsionDir = true;
} else {
return;
}
// if bond1 or bond2 are not to the stereo-controlling atoms, flip
// our expections of the torsion dir
int bond1AtomIdx = bond1->getOtherAtomIdx(dblBond->getBeginAtomIdx());
if (bond1AtomIdx != dblBond->getStereoAtoms()[0] &&
bond1AtomIdx != dblBond->getStereoAtoms()[1]) {
sameTorsionDir = !sameTorsionDir;
}
int bond2AtomIdx = bond2->getOtherAtomIdx(dblBond->getEndAtomIdx());
if (bond2AtomIdx != dblBond->getStereoAtoms()[0] &&
bond2AtomIdx != dblBond->getStereoAtoms()[1]) {
sameTorsionDir = !sameTorsionDir;
}
}
/*
Time for some clarificatory text, because this gets really
confusing really fast.
The dihedral angle analysis above is based on viewing things
with an atom order as follows:
1
\
2 = 3
\
4
so dihedrals > 90 correspond to sameDir=true
however, the stereochemistry representation is
based on something more like this:
2
\
1 = 3
\
4
(i.e. we consider the direction-setting single bonds to be
starting at the double-bonded atom)
*/
bool reverseBondDir = sameTorsionDir;
Atom *atom1 = dblBond->getBeginAtom(), *atom2 = dblBond->getEndAtom();
if (needsDir[bond1->getIdx()]) {
for (auto bidx : singleBondNbrs[bond1->getIdx()]) {
// std::cerr << " neighbor from: " << bond1->getIdx() << " " << bidx
// << ": " << needsDir[bidx] << std::endl;
if (needsDir[bidx]) {
followupBonds.push_back(mol.getBondWithIdx(bidx));
}
}
}
if (needsDir[bond2->getIdx()]) {
for (auto bidx : singleBondNbrs[bond2->getIdx()]) {
// std::cerr << " neighbor from: " << bond2->getIdx() << " " << bidx
// << ": " << needsDir[bidx] << std::endl;
if (needsDir[bidx]) {
followupBonds.push_back(mol.getBondWithIdx(bidx));
}
}
}
if (!needsDir[bond1->getIdx()]) {
if (!needsDir[bond2->getIdx()]) {
// check that we agree
} else {
if (bond1->getBeginAtom() != atom1) {
reverseBondDir = !reverseBondDir;
}
setBondDirRelativeToAtom(bond2, atom2, bond1->getBondDir(),
reverseBondDir, needsDir);
}
} else if (!needsDir[bond2->getIdx()]) {
if (bond2->getBeginAtom() != atom2) {
reverseBondDir = !reverseBondDir;
}
setBondDirRelativeToAtom(bond1, atom1, bond2->getBondDir(), reverseBondDir,
needsDir);
} else {
setBondDirRelativeToAtom(bond1, atom1, Bond::ENDDOWNRIGHT, false, needsDir);
setBondDirRelativeToAtom(bond2, atom2, Bond::ENDDOWNRIGHT, reverseBondDir,
needsDir);
}
needsDir[bond1->getIdx()] = 0;
needsDir[bond2->getIdx()] = 0;
if (obond1 && needsDir[obond1->getIdx()]) {
setBondDirRelativeToAtom(obond1, atom1, bond1->getBondDir(),
bond1->getBeginAtom() == atom1, needsDir);
needsDir[obond1->getIdx()] = 0;
}
if (obond2 && needsDir[obond2->getIdx()]) {
setBondDirRelativeToAtom(obond2, atom2, bond2->getBondDir(),
bond2->getBeginAtom() == atom2, needsDir);
needsDir[obond2->getIdx()] = 0;
}
#if 0
std::cerr << " 1:" << bond1->getIdx() << " ";
if (obond1)
std::cerr << obond1->getIdx() << std::endl;
else
std::cerr << "N/A" << std::endl;
std::cerr << " 2:" << bond2->getIdx() << " ";
if (obond2)
std::cerr << obond2->getIdx() << std::endl;
else
std::cerr << "N/A" << std::endl;
std::cerr << "**********************\n";
std::cerr << "**********************\n";
std::cerr << "**********************\n";
#endif
for (Bond *oDblBond : followupBonds) {
// std::cerr << "FOLLOWUP: " << oDblBond->getIdx() << " "
// << needsDir[oDblBond->getIdx()] << std::endl;
updateDoubleBondNeighbors(mol, oDblBond, conf, needsDir, singleBondCounts,
singleBondNbrs);
}
}
bool isBondCandidateForStereo(const Bond *bond) {
PRECONDITION(bond, "no bond");
if (bond->getBondType() == Bond::DOUBLE &&
bond->getStereo() != Bond::STEREOANY &&
bond->getBondDir() != Bond::EITHERDOUBLE &&
bond->getBeginAtom()->getDegree() > 1 &&
bond->getEndAtom()->getDegree() > 1 &&
shouldDetectDoubleBondStereo(bond)) {
return true;
}
return false;
}
const Atom *findHighestCIPNeighbor(const Atom *atom, const Atom *skipAtom) {
PRECONDITION(atom, "bad atom");
unsigned bestCipRank = 0;
const Atom *bestCipRankedAtom = nullptr;
const auto &mol = atom->getOwningMol();
for (const auto &index :
boost::make_iterator_range(mol.getAtomNeighbors(atom))) {
const auto neighbor = mol[index];
if (neighbor == skipAtom) {
continue;
}
unsigned cip = 0;
if (!neighbor->getPropIfPresent(common_properties::_CIPRank, cip)) {
// If at least one of the atoms doesn't have a CIP rank, the highest rank
// does not make sense, so return a nullptr.
return nullptr;
} else if (cip > bestCipRank || bestCipRankedAtom == nullptr) {
bestCipRank = cip;
bestCipRankedAtom = neighbor;
} else if (cip == bestCipRank) {
// This also doesn't make sense if there is a tie (if that's possible).
// We still keep the best CIP rank in case something better comes around
// (also not sure if that's possible).
BOOST_LOG(rdWarningLog)
<< "Warning: duplicate CIP ranks found in findHighestCIPNeighbor()"
<< std::endl;
bestCipRankedAtom = nullptr;
}
}
return bestCipRankedAtom;
}
} // namespace
namespace Chirality {
typedef std::pair<int, int> INT_PAIR;
typedef std::vector<INT_PAIR> INT_PAIR_VECT;
typedef std::vector<INT_PAIR>::iterator INT_PAIR_VECT_I;
typedef std::vector<INT_PAIR>::const_iterator INT_PAIR_VECT_CI;
typedef INT_VECT CIP_ENTRY;
typedef std::vector<CIP_ENTRY> CIP_ENTRY_VECT;
template <typename T>
void debugVect(const std::vector<T> arg) {
typename std::vector<T>::const_iterator viIt;
std::stringstream outS;
for (viIt = arg.begin(); viIt != arg.end(); viIt++) {
outS << *viIt << " ";
}
BOOST_LOG(rdDebugLog) << outS.str() << std::endl;
}
// --------------------------------------------------
//
// Calculates chiral invariants for the atoms of a molecule
// These are based on Labute's proposal in:
// "An Efficient Algorithm for the Determination of Topological
// RS Chirality" Journal of the CCG (1996)
//
// --------------------------------------------------
void buildCIPInvariants(const ROMol &mol, DOUBLE_VECT &res) {
PRECONDITION(res.size() >= mol.getNumAtoms(), "res vect too small");
int atsSoFar = 0;
//
// NOTE:
// If you make modifications to this, keep in mind that it is
// essential that the initial comparison of ranks behave properly.
// So, though it seems like it would makes sense to include
// information about the number of Hs (or charge, etc) in the CIP
// invariants, this will result in bad rankings. For example, in
// this molecule: OC[C@H](C)O, including the number of Hs would
// cause the methyl group (atom 3) to be ranked higher than the CH2
// connected to O (atom 1). This is totally wrong.
//
// We also don't include any pre-existing stereochemistry information.
// Though R and S assignments do factor in to the priorities of atoms,
// we're starting here from scratch and we'll let the R and S stuff
// be taken into account during the iterations.
//
for (ROMol::ConstAtomIterator atIt = mol.beginAtoms(); atIt != mol.endAtoms();
++atIt) {
const unsigned short nMassBits = 10;
const unsigned short maxMass = 1 << nMassBits;
Atom const *atom = *atIt;
unsigned long invariant = 0;
int num = atom->getAtomicNum() % 128;
// get an int with the deviation in the mass from the default:
int mass = 0;
if (atom->getIsotope()) {
mass =
atom->getIsotope() -
PeriodicTable::getTable()->getMostCommonIsotope(atom->getAtomicNum());
if (mass >= 0) {
mass += 1;
}
}
mass += maxMass / 2;
if (mass < 0) {
mass = 0;
} else {
mass = mass % maxMass;
}
#if 0
// NOTE: the inclusion of hybridization in the invariant (as
// suggested in the original paper), leads to the situation
// that
// C[C@@](O)(C=C)C(C)CC
// and
// C[C@@](O)(C=C)C(C)CO
// are assigned S chirality even though the rest of the world
// seems to agree that they ought to be R (atom 3, sp2, is ranked
// higher than atom 5, sp3, no matter what their environments)
int hyb=0;
switch(atom->getHybridization()) {
case Atom::SP: hyb=6;break;
case Atom::SP2: hyb=5;break;
case Atom::SP3: hyb=1;break;
case Atom::SP3D: hyb=3;break;
case Atom::SP3D2: hyb=2;break;
default: break;
}
#endif
invariant = num; // 7 bits here
invariant = (invariant << nMassBits) | mass;
int mapnum = -1;
atom->getPropIfPresent(common_properties::molAtomMapNumber, mapnum);
mapnum = (mapnum + 1) % 1024; // increment to allow map numbers of zero
// (though that would be stupid)
invariant = (invariant << 10) | mapnum;
res[atsSoFar++] = invariant;
}
}
void iterateCIPRanks(const ROMol &mol, const DOUBLE_VECT &invars,
UINT_VECT &ranks, bool seedWithInvars) {
PRECONDITION(invars.size() == mol.getNumAtoms(), "bad invars size");
PRECONDITION(ranks.size() >= mol.getNumAtoms(), "bad ranks size");
unsigned int numAtoms = mol.getNumAtoms();
CIP_ENTRY_VECT cipEntries(numAtoms);
for (auto &vec : cipEntries) {
vec.reserve(16);
}
#ifdef VERBOSE_CANON
BOOST_LOG(rdDebugLog) << "invariants:" << std::endl;
for (unsigned int i = 0; i < numAtoms; i++) {
BOOST_LOG(rdDebugLog) << i << ": " << invars[i] << std::endl;
}
#endif
// rank those:
Rankers::rankVect(invars, ranks);
#ifdef VERBOSE_CANON
BOOST_LOG(rdDebugLog) << "initial ranks:" << std::endl;
for (unsigned int i = 0; i < numAtoms; ++i) {
BOOST_LOG(rdDebugLog) << i << ": " << ranks[i] << std::endl;
}
#endif
// Start each atom's rank vector with its atomic number:
// Note: in general one should avoid the temptation to
// use invariants here, those lead to incorrect answers
for (unsigned int i = 0; i < numAtoms; i++) {
if (seedWithInvars) {
cipEntries[i].push_back(static_cast<int>(invars[i]));
} else {
cipEntries[i].push_back(mol[i]->getAtomicNum());
cipEntries[i].push_back(static_cast<int>(ranks[i]));
}
}
// Loop until either:
// 1) all classes are uniquified
// 2) the number of ranks doesn't change from one iteration to
// the next
// 3) we've gone through maxIts times
// maxIts is calculated by dividing the number of atoms
// by 2. That's a pessimal version of the
// maximum number of steps required for two atoms to
// "feel" each other (each influences one additional
// neighbor shell per iteration).
unsigned int maxIts = numAtoms / 2 + 1;
unsigned int numIts = 0;
int lastNumRanks = -1;
unsigned int numRanks = *std::max_element(ranks.begin(), ranks.end()) + 1;
std::vector<unsigned int> counts(ranks.size());
std::vector<unsigned int> updatedNbrIdxs;
updatedNbrIdxs.reserve(8);
while (numRanks < numAtoms && numIts < maxIts &&
(lastNumRanks < 0 ||
static_cast<unsigned int>(lastNumRanks) < numRanks)) {
unsigned int longestEntry = 0;
// ----------------------------------------------------
//
// for each atom, get a sorted list of its neighbors' ranks:
//
for (unsigned int index = 0; index < numAtoms; ++index) {
// Note: counts is cleaned up when we drain into cipEntries.
updatedNbrIdxs.clear();
// start by pushing on our neighbors' ranks:
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(mol[index]);
while (beg != end) {
const Bond *bond = mol[*beg];
++beg;
unsigned int nbrIdx = bond->getOtherAtomIdx(index);
updatedNbrIdxs.push_back(nbrIdx);
// put the neighbor in 2N times where N is the bond order as a double.
// this is to treat aromatic linkages on fair footing. i.e. at least in
// the first iteration --c(:c):c and --C(=C)-C should look the same.
// this was part of issue 3009911
// a special case for chiral phosphorus compounds
// (this was leading to incorrect assignment of R/S labels ):
bool isChiralPhosphorusSpecialCase = false;
if (bond->getBondType() == Bond::DOUBLE) {
const Atom *nbr = mol[nbrIdx];
if (nbr->getAtomicNum() == 15) {
unsigned int nbrDeg = nbr->getDegree();
isChiralPhosphorusSpecialCase = nbrDeg == 3 || nbrDeg == 4;
}
};
// general justification of this is:
// Paragraph 2.2. in the 1966 article is "Valence-Bond Conventions:
// Multiple-Bond Unsaturation and Aromaticity". It contains several
// conventions of which convention (b) is the one applying here:
// "(b) Contributions by d orbitals to bonds of quadriligant atoms are
// neglected."
// FIX: this applies to more than just P
if (isChiralPhosphorusSpecialCase) {
counts[nbrIdx] += 1;
} else {
counts[nbrIdx] += getTwiceBondType(*bond);
}
}
// For each of our neighbors' ranks weighted by bond type, copy it N times
// to our cipEntry in reverse rank order, where N is the weight.
if (updatedNbrIdxs.size() > 1) { // compare vs 1 for performance.
std::sort(std::begin(updatedNbrIdxs), std::end(updatedNbrIdxs),
[&ranks](unsigned int idx1, unsigned int idx2) {
return ranks[idx1] > ranks[idx2];
});
}
auto &cipEntry = cipEntries[index];
for (auto nbrIdx : updatedNbrIdxs) {
unsigned int count = counts[nbrIdx];
cipEntry.insert(cipEntry.end(), count, ranks[nbrIdx] + 1);
counts[nbrIdx] = 0;
}
// add a zero for each coordinated H as long as we're not a query atom
if (!mol[index]->hasQuery()) {
cipEntry.insert(cipEntry.end(), mol[index]->getTotalNumHs(), 0);
}
if (cipEntry.size() > longestEntry) {
longestEntry = rdcast<unsigned int>(cipEntry.size());
}
}
// ----------------------------------------------------
//
// pad the entries so that we compare rounds to themselves:
//
for (unsigned int index = 0; index < numAtoms; ++index) {
auto sz = rdcast<unsigned int>(cipEntries[index].size());
if (sz < longestEntry) {
cipEntries[index].insert(cipEntries[index].end(), longestEntry - sz,
-1);
}
}
// ----------------------------------------------------
//
// sort the new ranks and update the list of active indices:
//
lastNumRanks = numRanks;
Rankers::rankVect(cipEntries, ranks);
numRanks = *std::max_element(ranks.begin(), ranks.end()) + 1;
// now truncate each vector and stick the rank at the end
for (unsigned int i = 0; i < numAtoms; ++i) {
cipEntries[i][numIts + 1] = ranks[i];
cipEntries[i].erase(cipEntries[i].begin() + numIts + 2,
cipEntries[i].end());
}
++numIts;
#ifdef VERBOSE_CANON
BOOST_LOG(rdDebugLog) << "strings and ranks:" << std::endl;
for (unsigned int i = 0; i < numAtoms; i++) {
BOOST_LOG(rdDebugLog) << i << ": " << ranks[i] << " > ";
debugVect(cipEntries[i]);
}
#endif
}
}
// Figure out the CIP ranks for the atoms of a molecule
void assignAtomCIPRanks(const ROMol &mol, UINT_VECT &ranks) {
PRECONDITION((!ranks.size() || ranks.size() >= mol.getNumAtoms()),
"bad ranks size");
if (!ranks.size()) {
ranks.resize(mol.getNumAtoms());
}
unsigned int numAtoms = mol.getNumAtoms();
#ifndef USE_NEW_STEREOCHEMISTRY
// get the initial invariants:
DOUBLE_VECT invars(numAtoms, 0);
buildCIPInvariants(mol, invars);
iterateCIPRanks(mol, invars, ranks, false);
#else
Canon::chiralRankMolAtoms(mol, ranks);
#endif
// copy the ranks onto the atoms:
for (unsigned int i = 0; i < numAtoms; ++i) {
mol[i]->setProp(common_properties::_CIPRank, ranks[i], 1);
}
}
// construct a vector with <atomIdx,direction> pairs for
// neighbors of a given atom. This list will only be
// non-empty if at least one of the bonds has its direction
// set.
void findAtomNeighborDirHelper(const ROMol &mol, const Atom *atom,
const Bond *refBond, UINT_VECT &ranks,
INT_PAIR_VECT &neighbors,
bool &hasExplicitUnknownStereo) {
PRECONDITION(atom, "bad atom");
PRECONDITION(refBond, "bad bond");
bool seenDir = false;
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(atom);
while (beg != end) {
const Bond *bond = mol[*beg];
// check whether this bond is explicitly set to have unknown stereo
if (!hasExplicitUnknownStereo) {
int explicit_unknown_stereo;
if (bond->getBondDir() == Bond::UNKNOWN // there's a squiggle bond
|| (bond->getPropIfPresent<int>(common_properties::_UnknownStereo,
explicit_unknown_stereo) &&
explicit_unknown_stereo)) {
hasExplicitUnknownStereo = true;
}
}
Bond::BondDir dir = bond->getBondDir();
if (bond->getIdx() != refBond->getIdx()) {
if (dir == Bond::ENDDOWNRIGHT || dir == Bond::ENDUPRIGHT) {
seenDir = true;
// If we're considering the bond "backwards", (i.e. from end
// to beginning, reverse the effective direction:
if (atom != bond->getBeginAtom()) {
if (dir == Bond::ENDDOWNRIGHT) {
dir = Bond::ENDUPRIGHT;
} else {
dir = Bond::ENDDOWNRIGHT;
}
}
}
Atom *nbrAtom = bond->getOtherAtom(atom);
neighbors.push_back(std::make_pair(nbrAtom->getIdx(), dir));
}
++beg;
}
if (!seenDir) {
neighbors.clear();
} else {
if (neighbors.size() == 2 &&
ranks[neighbors[0].first] == ranks[neighbors[1].first]) {
// the two substituents are identical, no stereochemistry here:
neighbors.clear();
} else {
// it's possible that direction was set only one of the bonds, set the
// other
// bond's direction to be reversed:
if (neighbors[0].second != Bond::ENDDOWNRIGHT &&
neighbors[0].second != Bond::ENDUPRIGHT) {
CHECK_INVARIANT(neighbors.size() > 1, "too few neighbors");
neighbors[0].second = neighbors[1].second == Bond::ENDDOWNRIGHT
? Bond::ENDUPRIGHT
: Bond::ENDDOWNRIGHT;
} else if (neighbors.size() > 1 &&
neighbors[1].second != Bond::ENDDOWNRIGHT &&
neighbors[1].second != Bond::ENDUPRIGHT) {
neighbors[1].second = neighbors[0].second == Bond::ENDDOWNRIGHT
? Bond::ENDUPRIGHT
: Bond::ENDDOWNRIGHT;
}
}
}
}
// find the neighbors for an atoms that are not connected by single bond that is
// not refBond
// if checkDir is true only neighbor atoms with bonds marked with a direction
// will be returned
void findAtomNeighborsHelper(const ROMol &mol, const Atom *atom,
const Bond *refBond, UINT_VECT &neighbors,
bool checkDir = false,
bool includeAromatic = false) {
PRECONDITION(atom, "bad atom");
PRECONDITION(refBond, "bad bond");
neighbors.clear();
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(atom);
while (beg != end) {
const Bond *bond = mol[*beg];
Bond::BondDir dir = bond->getBondDir();
if ((bond->getBondType() == Bond::SINGLE ||
(includeAromatic && bond->getBondType() == Bond::AROMATIC)) &&
bond->getIdx() != refBond->getIdx()) {
if (checkDir) {
if ((dir != Bond::ENDDOWNRIGHT) && (dir != Bond::ENDUPRIGHT)) {
++beg;
continue;
}
}
Atom *nbrAtom = bond->getOtherAtom(atom);
neighbors.push_back(nbrAtom->getIdx());
}
++beg;
}
}
// conditions for an atom to be a candidate for ring stereochem:
// 1) two non-ring neighbors that have different ranks
// 2) one non-ring neighbor and two ring neighbors (the ring neighbors will
// have the same rank)
// 3) four ring neighbors with three different ranks
// 4) three ring neighbors with two different ranks
// example for this last one: C[C@H]1CC2CCCC3CCCC(C1)[C@@H]23
// Note that N atoms are only candidates if they are in a 3-ring
bool atomIsCandidateForRingStereochem(const ROMol &mol, const Atom *atom) {
PRECONDITION(atom, "bad atom");
bool res = false;
std::set<unsigned int> nbrRanks;
if (!atom->getPropIfPresent(common_properties::_ringStereochemCand, res)) {
const RingInfo *ringInfo = mol.getRingInfo();
if (ringInfo->isInitialized() && ringInfo->numAtomRings(atom->getIdx())) {
// three-coordinate N additional requirements:
// in a ring of size 3 (from InChI)
// OR
// a bridgehead (RDKit extension)
if (atom->getAtomicNum() == 7 && atom->getDegree() == 3 &&
!ringInfo->isAtomInRingOfSize(atom->getIdx(), 3) &&
!queryIsAtomBridgehead(atom)) {
return false;
}
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(atom);
std::vector<const Atom *> nonRingNbrs;
std::vector<const Atom *> ringNbrs;
while (beg != end) {
const Bond *bond = mol[*beg];
if (!ringInfo->numBondRings(bond->getIdx())) {
nonRingNbrs.push_back(bond->getOtherAtom(atom));
} else {
const Atom *nbr = bond->getOtherAtom(atom);
ringNbrs.push_back(nbr);
unsigned int rnk = 0;
nbr->getPropIfPresent(common_properties::_CIPRank, rnk);
nbrRanks.insert(rnk);
}
++beg;
}
unsigned int rank1 = 0, rank2 = 0;
switch (nonRingNbrs.size()) {
case 2:
if (nonRingNbrs[0]->getPropIfPresent(common_properties::_CIPRank,
rank1) &&
nonRingNbrs[1]->getPropIfPresent(common_properties::_CIPRank,
rank2)) {
if (rank1 == rank2) {
res = false;
} else {
res = true;
}
}
break;
case 1:
if (ringNbrs.size() >= 2) {
res = true;
}
break;
case 0:
if (ringNbrs.size() == 4 && nbrRanks.size() == 3) {
res = true;
} else if (ringNbrs.size() == 3 && nbrRanks.size() == 2) {
res = true;
} else {
res = false;
}
break;
default:
res = false;
}
}
atom->setProp(common_properties::_ringStereochemCand, res, 1);
}
return res;
}
// finds all possible chiral special cases.
// at the moment this is just candidates for ring stereochemistry
void findChiralAtomSpecialCases(ROMol &mol,
boost::dynamic_bitset<> &possibleSpecialCases) {
PRECONDITION(possibleSpecialCases.size() >= mol.getNumAtoms(),
"bit vector too small");
possibleSpecialCases.reset();
if (!mol.getRingInfo()->isInitialized()) {
VECT_INT_VECT sssrs;
MolOps::symmetrizeSSSR(mol, sssrs);
}
boost::dynamic_bitset<> atomsSeen(mol.getNumAtoms());
boost::dynamic_bitset<> atomsUsed(mol.getNumAtoms());
boost::dynamic_bitset<> bondsSeen(mol.getNumBonds());
for (ROMol::AtomIterator ait = mol.beginAtoms(); ait != mol.endAtoms();
++ait) {
const Atom *atom = *ait;
if (atomsSeen[atom->getIdx()]) {
continue;
}
if (atom->getChiralTag() == Atom::CHI_UNSPECIFIED ||
atom->hasProp(common_properties::_CIPCode) ||
!mol.getRingInfo()->numAtomRings(atom->getIdx()) ||
!atomIsCandidateForRingStereochem(mol, atom)) {
continue;
}
// do a BFS from this ring atom along ring bonds and find other
// stereochemistry candidates.
std::list<const Atom *> nextAtoms;
// start with finding viable neighbors
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(atom);
while (beg != end) {
unsigned int bidx = mol[*beg]->getIdx();
if (!bondsSeen[bidx]) {
bondsSeen.set(bidx);
if (mol.getRingInfo()->numBondRings(bidx)) {
const Atom *oatom = mol[*beg]->getOtherAtom(atom);
if (!atomsSeen[oatom->getIdx()]) {
nextAtoms.push_back(oatom);
atomsUsed.set(oatom->getIdx());
}
}
}
++beg;
}
INT_VECT ringStereoAtoms(0);
if (!nextAtoms.empty()) {
atom->getPropIfPresent(common_properties::_ringStereoAtoms,
ringStereoAtoms);
}
while (!nextAtoms.empty()) {
const Atom *ratom = nextAtoms.front();
nextAtoms.pop_front();
atomsSeen.set(ratom->getIdx());
if (ratom->getChiralTag() != Atom::CHI_UNSPECIFIED &&
!ratom->hasProp(common_properties::_CIPCode) &&
atomIsCandidateForRingStereochem(mol, ratom)) {
int same = (ratom->getChiralTag() == atom->getChiralTag()) ? 1 : -1;
ringStereoAtoms.push_back(same * (ratom->getIdx() + 1));
INT_VECT oringatoms(0);
ratom->getPropIfPresent(common_properties::_ringStereoAtoms,
oringatoms);
oringatoms.push_back(same * (atom->getIdx() + 1));
ratom->setProp(common_properties::_ringStereoAtoms, oringatoms, true);
possibleSpecialCases.set(ratom->getIdx());
possibleSpecialCases.set(atom->getIdx());
}
// now push this atom's neighbors
boost::tie(beg, end) = mol.getAtomBonds(ratom);
while (beg != end) {
unsigned int bidx = mol[*beg]->getIdx();
if (!bondsSeen[bidx]) {
bondsSeen.set(bidx);
if (mol.getRingInfo()->numBondRings(bidx)) {
const Atom *oatom = mol[*beg]->getOtherAtom(ratom);
if (!atomsSeen[oatom->getIdx()] && !atomsUsed[oatom->getIdx()]) {
nextAtoms.push_back(oatom);
atomsUsed.set(oatom->getIdx());
}
}
}
++beg;
}
} // end of BFS
if (ringStereoAtoms.size() != 0) {
atom->setProp(common_properties::_ringStereoAtoms, ringStereoAtoms, true);
// because we're only going to hit each ring atom once, the first atom we
// encounter in a ring is going to end up with all the other atoms set as
// stereoAtoms, but each of them will only have the first atom present. We
// need to fix that. because the traverse from the first atom only
// followed ring bonds, these things are all by definition in one ring
// system. (Q: is this true if there's a spiro center in there?)
INT_VECT same(mol.getNumAtoms(), 0);
for (auto ringAtomEntry : ringStereoAtoms) {
int ringAtomIdx =
ringAtomEntry < 0 ? -ringAtomEntry - 1 : ringAtomEntry - 1;
same[ringAtomIdx] = ringAtomEntry;
}
for (INT_VECT_CI rae = ringStereoAtoms.begin();
rae != ringStereoAtoms.end(); ++rae) {
int ringAtomEntry = *rae;
int ringAtomIdx =
ringAtomEntry < 0 ? -ringAtomEntry - 1 : ringAtomEntry - 1;
INT_VECT lringatoms(0);
mol.getAtomWithIdx(ringAtomIdx)
->getPropIfPresent(common_properties::_ringStereoAtoms, lringatoms);
CHECK_INVARIANT(lringatoms.size() > 0, "no other ring atoms found.");
for (auto orae = rae + 1; orae != ringStereoAtoms.end(); ++orae) {
int oringAtomEntry = *orae;
int oringAtomIdx =
oringAtomEntry < 0 ? -oringAtomEntry - 1 : oringAtomEntry - 1;
int theseDifferent = (ringAtomEntry < 0) ^ (oringAtomEntry < 0);
lringatoms.push_back(theseDifferent ? -(oringAtomIdx + 1)
: (oringAtomIdx + 1));
INT_VECT olringatoms(0);
mol.getAtomWithIdx(oringAtomIdx)
->getPropIfPresent(common_properties::_ringStereoAtoms,
olringatoms);
CHECK_INVARIANT(olringatoms.size() > 0, "no other ring atoms found.");
olringatoms.push_back(theseDifferent ? -(ringAtomIdx + 1)
: (ringAtomIdx + 1));
mol.getAtomWithIdx(oringAtomIdx)
->setProp(common_properties::_ringStereoAtoms, olringatoms);
}
mol.getAtomWithIdx(ringAtomIdx)
->setProp(common_properties::_ringStereoAtoms, lringatoms);
}
} else {
possibleSpecialCases.reset(atom->getIdx());
}
atomsSeen.set(atom->getIdx());
}
}
std::pair<bool, bool> isAtomPotentialChiralCenter(
const Atom *atom, const ROMol &mol, const UINT_VECT &ranks,
Chirality::INT_PAIR_VECT &nbrs) {
// loop over all neighbors and form a decorated list of their
// ranks:
bool legalCenter = true;
bool hasDupes = false;
if (atom->getTotalDegree() > 4) {
// we only know tetrahedral chirality
legalCenter = false;
} else {
// cases we can exclude immediately without having to look at neighbors
// ranks:
if (atom->getTotalDegree() < 3) {
legalCenter = false;
} else if (atom->getDegree() < 3 &&
(atom->getAtomicNum() != 15 && atom->getAtomicNum() != 33)) {
// less than three neighbors is never stereogenic
// unless it is a phosphine/arsine with implicit H (this is from InChI)
legalCenter = false;
} else if (atom->getDegree() == 3 && atom->getTotalNumHs() != 1) {
// assume something that's really three coordinate isn't potentially
// chiral, then look for exceptions
legalCenter = false;
if (atom->getAtomicNum() == 7) {
// three-coordinate N additional requirements:
// in a ring of size 3 (from InChI)
// OR
/// is a bridgehead atom (RDKit extension)
if (mol.getRingInfo()->isAtomInRingOfSize(atom->getIdx(), 3) ||
queryIsAtomBridgehead(atom)) {
legalCenter = true;
}
} else if (atom->getAtomicNum() == 15 || atom->getAtomicNum() == 33) {
// three-coordinate phosphines and arsines
// are always treated as stereogenic even with H atom neighbors.
// (this is from InChI)
legalCenter = true;
} else if (atom->getAtomicNum() == 16 || atom->getAtomicNum() == 34) {
if (atom->getExplicitValence() == 4 ||
(atom->getExplicitValence() == 3 && atom->getFormalCharge() == 1)) {
// we also accept sulfur or selenium with either a positive charge
// or a double bond:
legalCenter = true;
}
}
}
if (legalCenter) {
boost::dynamic_bitset<> codesSeen(mol.getNumAtoms());
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(atom);
while (beg != end) {
unsigned int otherIdx = mol[*beg]->getOtherAtom(atom)->getIdx();
CHECK_INVARIANT(ranks[otherIdx] < mol.getNumAtoms(),
"CIP rank higher than the number of atoms.");
// watch for neighbors with duplicate ranks, which would mean
// that we cannot be chiral:
if (codesSeen[ranks[otherIdx]]) {
// we've already seen this code, it's a dupe
hasDupes = true;
break;
}
codesSeen[ranks[otherIdx]] = 1;
nbrs.push_back(std::make_pair(ranks[otherIdx], mol[*beg]->getIdx()));
++beg;
}
}
}
return std::make_pair(legalCenter, hasDupes);
}
// returns a pair:
// 1) are there unassigned stereoatoms
// 2) did we assign any?
std::pair<bool, bool> assignAtomChiralCodes(ROMol &mol, UINT_VECT &ranks,
bool flagPossibleStereoCenters) {
PRECONDITION((!ranks.size() || ranks.size() == mol.getNumAtoms()),
"bad rank vector size");
bool atomChanged = false;
unsigned int unassignedAtoms = 0;
// ------------------
// now loop over each atom and, if it's marked as chiral,
// figure out the appropriate CIP label:
for (ROMol::AtomIterator atIt = mol.beginAtoms(); atIt != mol.endAtoms();
++atIt) {
Atom *atom = *atIt;
Atom::ChiralType tag = atom->getChiralTag();
// only worry about this atom if it has a marked chirality
// we understand:
if (flagPossibleStereoCenters ||
(tag != Atom::CHI_UNSPECIFIED && tag != Atom::CHI_OTHER)) {
if (atom->hasProp(common_properties::_CIPCode)) {
continue;
}
if (!ranks.size()) {
// if we need to, get the "CIP" ranking of each atom:
assignAtomCIPRanks(mol, ranks);
}
Chirality::INT_PAIR_VECT nbrs;
bool legalCenter, hasDupes;
// note that hasDupes is only evaluated if legalCenter==true
boost::tie(legalCenter, hasDupes) =
isAtomPotentialChiralCenter(atom, mol, ranks, nbrs);
if (legalCenter) {
++unassignedAtoms;
}
if (legalCenter && !hasDupes && flagPossibleStereoCenters) {
atom->setProp(common_properties::_ChiralityPossible, 1);
}
if (legalCenter && !hasDupes && tag != Atom::CHI_UNSPECIFIED &&
tag != Atom::CHI_OTHER) {
// stereochem is possible and we have no duplicate neighbors, assign
// a CIP code:
atomChanged = true;
--unassignedAtoms;
// sort the list of neighbors by their CIP ranks:
std::sort(nbrs.begin(), nbrs.end(), Rankers::pairLess<int, int>());
// collect the list of neighbor indices:
std::list<int> nbrIndices;
for (Chirality::INT_PAIR_VECT_CI nbrIt = nbrs.begin();
nbrIt != nbrs.end(); ++nbrIt) {
nbrIndices.push_back((*nbrIt).second);
}
// ask the atom how many swaps we have to make:
int nSwaps = atom->getPerturbationOrder(nbrIndices);
// if the atom has 3 neighbors and a hydrogen, add a swap:
if (nbrIndices.size() == 3 && atom->getTotalNumHs() == 1) {
++nSwaps;
}
// if that number is odd, we'll change our chirality:
if (nSwaps % 2) {
if (tag == Atom::CHI_TETRAHEDRAL_CCW) {
tag = Atom::CHI_TETRAHEDRAL_CW;
} else {
tag = Atom::CHI_TETRAHEDRAL_CCW;
}
}
// now assign the CIP code:
std::string cipCode;
if (tag == Atom::CHI_TETRAHEDRAL_CCW) {
cipCode = "S";
} else {
cipCode = "R";
}
atom->setProp(common_properties::_CIPCode, cipCode);
}
}
}
return std::make_pair((unassignedAtoms > 0), atomChanged);
}
// returns a pair:
// 1) are there unassigned stereo bonds?
// 2) did we assign any?
std::pair<bool, bool> assignBondStereoCodes(ROMol &mol, UINT_VECT &ranks) {
PRECONDITION((!ranks.size() || ranks.size() == mol.getNumAtoms()),
"bad rank vector size");
bool assignedABond = false;
unsigned int unassignedBonds = 0;
boost::dynamic_bitset<> bondsToClear(mol.getNumBonds());
// find the double bonds:
for (ROMol::BondIterator bondIt = mol.beginBonds(); bondIt != mol.endBonds();
++bondIt) {
if ((*bondIt)->getBondType() == Bond::DOUBLE) {
Bond *dblBond = *bondIt;
if (dblBond->getStereo() != Bond::STEREONONE) {
continue;
}
if (!ranks.size()) {
assignAtomCIPRanks(mol, ranks);
}
dblBond->getStereoAtoms().clear();
// at the moment we are ignoring stereochem on ring bonds with less than
// 8
// members.
if (shouldDetectDoubleBondStereo(dblBond)) {
const Atom *begAtom = dblBond->getBeginAtom();
const Atom *endAtom = dblBond->getEndAtom();
// we're only going to handle 2 or three coordinate atoms:
if ((begAtom->getDegree() == 2 || begAtom->getDegree() == 3) &&
(endAtom->getDegree() == 2 || endAtom->getDegree() == 3)) {
++unassignedBonds;
// look around each atom and see if it has at least one bond with
// direction marked:
// the pairs here are: atomrank,bonddir
Chirality::INT_PAIR_VECT begAtomNeighbors, endAtomNeighbors;
bool hasExplicitUnknownStereo = false;
int bgn_stereo = false, end_stereo = false;
if ((dblBond->getBeginAtom()->getPropIfPresent(
common_properties::_UnknownStereo, bgn_stereo) &&
bgn_stereo) ||
(dblBond->getEndAtom()->getPropIfPresent(
common_properties::_UnknownStereo, end_stereo) &&
end_stereo)) {
hasExplicitUnknownStereo = true;
}
Chirality::findAtomNeighborDirHelper(mol, begAtom, dblBond, ranks,
begAtomNeighbors,
hasExplicitUnknownStereo);
Chirality::findAtomNeighborDirHelper(mol, endAtom, dblBond, ranks,
endAtomNeighbors,
hasExplicitUnknownStereo);
if (begAtomNeighbors.size() && endAtomNeighbors.size()) {
// Each atom has at least one neighboring bond with marked
// directionality. Find the highest-ranked directionality
// on each side:
int begDir, endDir, endNbrAid, begNbrAid;
if (begAtomNeighbors.size() == 1 ||
ranks[begAtomNeighbors[0].first] >
ranks[begAtomNeighbors[1].first]) {
begDir = begAtomNeighbors[0].second;
begNbrAid = begAtomNeighbors[0].first;
} else {
begDir = begAtomNeighbors[1].second;
begNbrAid = begAtomNeighbors[1].first;
}
if (endAtomNeighbors.size() == 1 ||
ranks[endAtomNeighbors[0].first] >
ranks[endAtomNeighbors[1].first]) {
endDir = endAtomNeighbors[0].second;
endNbrAid = endAtomNeighbors[0].first;
} else {
endDir = endAtomNeighbors[1].second;
endNbrAid = endAtomNeighbors[1].first;
}
bool conflictingBegin =
(begAtomNeighbors.size() == 2 &&
begAtomNeighbors[0].second == begAtomNeighbors[1].second);
bool conflictingEnd =
(endAtomNeighbors.size() == 2 &&
endAtomNeighbors[0].second == endAtomNeighbors[1].second);
if (conflictingBegin || conflictingEnd) {
dblBond->setStereo(Bond::STEREONONE);
BOOST_LOG(rdWarningLog) << "Conflicting single bond directions "
"around double bond at index "
<< dblBond->getIdx() << "." << std::endl;
BOOST_LOG(rdWarningLog) << " BondStereo set to STEREONONE and "
"single bond directions set to NONE."
<< std::endl;
assignedABond = true;
if (conflictingBegin) {
bondsToClear[mol.getBondBetweenAtoms(begAtomNeighbors[0].first,
begAtom->getIdx())
->getIdx()] = 1;
bondsToClear[mol.getBondBetweenAtoms(begAtomNeighbors[1].first,
begAtom->getIdx())
->getIdx()] = 1;
}
if (conflictingEnd) {
bondsToClear[mol.getBondBetweenAtoms(endAtomNeighbors[0].first,
endAtom->getIdx())
->getIdx()] = 1;
bondsToClear[mol.getBondBetweenAtoms(endAtomNeighbors[1].first,
endAtom->getIdx())
->getIdx()] = 1;
}
} else {
dblBond->getStereoAtoms().push_back(begNbrAid);
dblBond->getStereoAtoms().push_back(endNbrAid);
if (hasExplicitUnknownStereo) {
dblBond->setStereo(Bond::STEREOANY);
assignedABond = true;
} else if (begDir == endDir) {
// In findAtomNeighborDirHelper, we've set up the
// bond directions here so that they correspond to
// having both single bonds START at the double bond.
// This means that if the single bonds point in the same
// direction, the bond is cis, "Z"
dblBond->setStereo(Bond::STEREOZ);
assignedABond = true;
} else {
dblBond->setStereo(Bond::STEREOE);
assignedABond = true;
}
}
--unassignedBonds;
}
}
}
}
}
for (unsigned int i = 0; i < mol.getNumBonds(); ++i) {
if (bondsToClear[i]) {
mol.getBondWithIdx(i)->setBondDir(Bond::NONE);
}
}
return std::make_pair(unassignedBonds > 0, assignedABond);
}
// reassign atom ranks by supplementing the current ranks
// with information about known chirality
void rerankAtoms(const ROMol &mol, UINT_VECT &ranks) {
PRECONDITION(ranks.size() == mol.getNumAtoms(), "bad rank vector size");
unsigned int factor = 100;
while (factor < mol.getNumAtoms()) {
factor *= 10;
}
#ifdef VERBOSE_CANON
BOOST_LOG(rdDebugLog) << "rerank PRE: " << std::endl;
for (int i = 0; i < mol.getNumAtoms(); i++) {
BOOST_LOG(rdDebugLog) << " " << i << ": " << ranks[i] << std::endl;
}
#endif
DOUBLE_VECT invars(mol.getNumAtoms());
// and now supplement them:
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
invars[i] = ranks[i] * factor;
const Atom *atom = mol.getAtomWithIdx(i);
// Priority order: R > S > nothing
std::string cipCode;
if (atom->getPropIfPresent(common_properties::_CIPCode, cipCode)) {
if (cipCode == "S") {
invars[i] += 10;
} else if (cipCode == "R") {
invars[i] += 20;
}
}
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(atom);
while (beg != end) {
const Bond *oBond = mol[*beg];
if (oBond->getBondType() == Bond::DOUBLE) {
if (oBond->getStereo() == Bond::STEREOE) {
invars[i] += 1;
} else if (oBond->getStereo() == Bond::STEREOZ) {
invars[i] += 2;
}
}
++beg;
}
}
iterateCIPRanks(mol, invars, ranks, true);
// copy the ranks onto the atoms:
for (unsigned int i = 0; i < mol.getNumAtoms(); i++) {
mol.getAtomWithIdx(i)->setProp(common_properties::_CIPRank, ranks[i]);
}
#ifdef VERBOSE_CANON
BOOST_LOG(rdDebugLog) << " post: " << std::endl;
for (int i = 0; i < mol.getNumAtoms(); i++) {
BOOST_LOG(rdDebugLog) << " " << i << ": " << ranks[i] << std::endl;
}
#endif
}
bool hasStereoBondDir(const Bond *bond) {
PRECONDITION(bond, "no bond");
return bond->getBondDir() == Bond::BondDir::ENDDOWNRIGHT ||
bond->getBondDir() == Bond::BondDir::ENDUPRIGHT;
}
const Bond *getNeighboringDirectedBond(const ROMol &mol, const Atom *atom) {
PRECONDITION(atom, "no atom");
for (const auto &bondIdx :
boost::make_iterator_range(mol.getAtomBonds(atom))) {
const Bond *bond = mol[bondIdx];
if (bond->getBondType() != Bond::BondType::DOUBLE &&
hasStereoBondDir(bond)) {
return bond;
}
}
return nullptr;
}
Bond::BondStereo translateEZLabelToCisTrans(Bond::BondStereo label) {
switch (label) {
case Bond::STEREOE:
return Bond::STEREOTRANS;
case Bond::STEREOZ:
return Bond::STEREOCIS;
default:
return label;
}
}
INT_VECT findStereoAtoms(const Bond *bond) {
PRECONDITION(bond, "bad bond");
PRECONDITION(bond->hasOwningMol(), "no mol");
PRECONDITION(bond->getBondType() == Bond::DOUBLE, "not double bond");
PRECONDITION(bond->getStereo() > Bond::BondStereo::STEREOANY,
"no defined stereo");
if (!bond->getStereoAtoms().empty()) {
return bond->getStereoAtoms();
}
if (bond->getStereo() == Bond::BondStereo::STEREOE ||
bond->getStereo() == Bond::BondStereo::STEREOZ) {
const Atom *startStereoAtom =
findHighestCIPNeighbor(bond->getBeginAtom(), bond->getEndAtom());
const Atom *endStereoAtom =
findHighestCIPNeighbor(bond->getEndAtom(), bond->getBeginAtom());
if (startStereoAtom == nullptr || endStereoAtom == nullptr) {
return {};
}
int startStereoAtomIdx = static_cast<int>(startStereoAtom->getIdx());
int endStereoAtomIdx = static_cast<int>(endStereoAtom->getIdx());
return {startStereoAtomIdx, endStereoAtomIdx};
} else {
BOOST_LOG(rdWarningLog) << "Unable to assign stereo atoms for bond "
<< bond->getIdx() << std::endl;
return {};
}
}
} // namespace Chirality
namespace MolOps {
/*
We're going to do this iteratively:
1) assign atom stereochemistry
2) assign bond stereochemistry
3) if there are still unresolved atoms or bonds
repeat the above steps as necessary
*/
void assignStereochemistry(ROMol &mol, bool cleanIt, bool force,
bool flagPossibleStereoCenters) {
if (!force && mol.hasProp(common_properties::_StereochemDone)) {
return;
}
// later we're going to need ring information, get it now if we don't
// have it already:
if (!mol.getRingInfo()->isInitialized()) {
MolOps::fastFindRings(mol);
}
#if 0
std::cerr << ">>>>>>>>>>>>>\n";
std::cerr << "assign stereochem\n";
mol.debugMol(std::cerr);
#endif
// as part of the preparation, we'll loop over the atoms and
// bonds to see if anything has stereochemistry
// indicated. There's no point in doing the work here if there
// are neither stereocenters nor bonds that we need to consider.
// The exception to this is when flagPossibleStereoCenters is
// true; then we always need to do the work
bool hasStereoAtoms = flagPossibleStereoCenters;
for (ROMol::AtomIterator atIt = mol.beginAtoms(); atIt != mol.endAtoms();
++atIt) {
if (cleanIt) {
if ((*atIt)->hasProp(common_properties::_CIPCode)) {
(*atIt)->clearProp(common_properties::_CIPCode);
}
if ((*atIt)->hasProp(common_properties::_ChiralityPossible)) {
(*atIt)->clearProp(common_properties::_ChiralityPossible);
}
}
if (!hasStereoAtoms && (*atIt)->getChiralTag() != Atom::CHI_UNSPECIFIED &&
(*atIt)->getChiralTag() != Atom::CHI_OTHER) {
hasStereoAtoms = true;
}
}
bool hasStereoBonds = false;
for (ROMol::BondIterator bondIt = mol.beginBonds(); bondIt != mol.endBonds();
++bondIt) {
if (cleanIt) {
if ((*bondIt)->getBondType() == Bond::DOUBLE) {
if ((*bondIt)->getBondDir() == Bond::EITHERDOUBLE) {
(*bondIt)->setStereo(Bond::STEREOANY);
} else if ((*bondIt)->getStereo() != Bond::STEREOANY) {
(*bondIt)->setStereo(Bond::STEREONONE);
(*bondIt)->getStereoAtoms().clear();
}
}
}
if (!hasStereoBonds && (*bondIt)->getBondType() == Bond::DOUBLE) {
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds((*bondIt)->getBeginAtom());
while (!hasStereoBonds && beg != end) {
const Bond *nbond = mol[*beg];
++beg;
if (nbond->getBondDir() == Bond::ENDDOWNRIGHT ||
nbond->getBondDir() == Bond::ENDUPRIGHT) {
hasStereoBonds = true;
}
}
boost::tie(beg, end) = mol.getAtomBonds((*bondIt)->getEndAtom());
while (!hasStereoBonds && beg != end) {
const Bond *nbond = mol[*beg];
++beg;
if (nbond->getBondDir() == Bond::ENDDOWNRIGHT ||
nbond->getBondDir() == Bond::ENDUPRIGHT) {
hasStereoBonds = true;
}
}
}
if (!cleanIt && hasStereoBonds) {
break; // no reason to keep iterating if we've already
// determined there are stereo bonds to consider
}
}
UINT_VECT atomRanks;
bool keepGoing = hasStereoAtoms | hasStereoBonds;
bool changedStereoAtoms, changedStereoBonds;
while (keepGoing) {
if (hasStereoAtoms) {
boost::tie(hasStereoAtoms, changedStereoAtoms) =
Chirality::assignAtomChiralCodes(mol, atomRanks,
flagPossibleStereoCenters);
} else {
changedStereoAtoms = false;
}
if (hasStereoBonds) {
boost::tie(hasStereoBonds, changedStereoBonds) =
Chirality::assignBondStereoCodes(mol, atomRanks);
} else {
changedStereoBonds = false;
}
keepGoing = (hasStereoAtoms || hasStereoBonds) &&
(changedStereoAtoms || changedStereoBonds);
if (keepGoing) {
// update the atom ranks based on the new information we have:
Chirality::rerankAtoms(mol, atomRanks);
}
#if 0
std::cout << "*************** done iteration " << keepGoing
<< " ***********" << std::endl;
mol.debugMol(std::cout);
std::cout << "*************** done iteration " << keepGoing
<< " ***********" << std::endl;
#endif
}
if (cleanIt) {
// if the ranks are needed again, this will force them to be
// re-calculated based on the stereo calculated above.
// atomRanks.clear();
for (ROMol::AtomIterator atIt = mol.beginAtoms(); atIt != mol.endAtoms();
++atIt) {
if ((*atIt)->hasProp(common_properties::_ringStereochemCand)) {
(*atIt)->clearProp(common_properties::_ringStereochemCand);
}
if ((*atIt)->hasProp(common_properties::_ringStereoAtoms)) {
(*atIt)->clearProp(common_properties::_ringStereoAtoms);
}
}
boost::dynamic_bitset<> possibleSpecialCases(mol.getNumAtoms());
Chirality::findChiralAtomSpecialCases(mol, possibleSpecialCases);
for (auto atom : mol.atoms()) {
if (atom->getChiralTag() != Atom::CHI_UNSPECIFIED &&
!atom->hasProp(common_properties::_CIPCode) &&
(!possibleSpecialCases[atom->getIdx()] ||
!atom->hasProp(common_properties::_ringStereoAtoms))) {
atom->setChiralTag(Atom::CHI_UNSPECIFIED);
// If the atom has an explicit hydrogen and no charge, that H
// was probably put there solely because of the chirality.
// So we'll go ahead and remove it.
// This was Issue 194
if (atom->getNumExplicitHs() == 1 && atom->getFormalCharge() == 0 &&
!atom->getIsAromatic()) {
atom->setNumExplicitHs(0);
atom->setNoImplicit(false);
atom->calcExplicitValence(false);
atom->calcImplicitValence(false);
}
}
}
for (auto bond : mol.bonds()) {
// wedged bonds to atoms that have no stereochem
// should be removed. (github issue 87)
if ((bond->getBondDir() == Bond::BEGINWEDGE ||
bond->getBondDir() == Bond::BEGINDASH) &&
bond->getBeginAtom()->getChiralTag() == Atom::CHI_UNSPECIFIED &&
bond->getEndAtom()->getChiralTag() == Atom::CHI_UNSPECIFIED) {
bond->setBondDir(Bond::NONE);
}
// check for directionality on single bonds around
// double bonds without stereo. This was github #2422
if (bond->getBondType() == Bond::DOUBLE &&
(bond->getStereo() == Bond::STEREOANY ||
bond->getStereo() == Bond::STEREONONE)) {
std::vector<Atom *> batoms = {bond->getBeginAtom(), bond->getEndAtom()};
for (auto batom : batoms) {
for (const auto &nbri :
boost::make_iterator_range(mol.getAtomBonds(batom))) {
auto nbrBndI = mol[nbri];
if ((nbrBndI->getBondDir() == Bond::ENDDOWNRIGHT ||
nbrBndI->getBondDir() == Bond::ENDUPRIGHT) &&
(nbrBndI->getBondType() == Bond::SINGLE ||
nbrBndI->getBondType() == Bond::AROMATIC)) {
// direction is set, and we know it's not because of our
// bond. What about other neighbors?
bool okToClear = true;
for (const auto &nbrj : boost::make_iterator_range(
mol.getAtomBonds(nbrBndI->getOtherAtom(batom)))) {
auto nbrBndJ = mol[nbrj];
if (nbrBndJ->getBondType() == Bond::DOUBLE &&
nbrBndJ->getStereo() != Bond::STEREOANY &&
nbrBndJ->getStereo() != Bond::STEREONONE) {
okToClear = false;
break;
}
}
if (okToClear) {
nbrBndI->setBondDir(Bond::NONE);
}
}
}
}
}
#if 0
// make sure CIS/TRANS assignments are actually stereo bonds
if ((*bondIt)->getBondType() == Bond::DOUBLE) {
if ((*bondIt)->getStereo() == Bond::STEREOCIS ||
(*bondIt)->getStereo() == Bond::STEREOTRANS) {
if (!atomRanks.size()) {
Chirality::assignAtomCIPRanks(mol, atomRanks);
}
const Atom *begAtom = (*bondIt)->getBeginAtom(),
*endAtom = (*bondIt)->getEndAtom();
UINT_VECT begAtomNeighbors, endAtomNeighbors;
Chirality::findAtomNeighborsHelper(mol, begAtom, *bondIt,
begAtomNeighbors);
Chirality::findAtomNeighborsHelper(mol, endAtom, *bondIt,
endAtomNeighbors);
// Note, this relies on this being a hydrogen-suppressed
// graph as the 'Note' in the doc string of this function
// indicates is a pre-condition.
if ((begAtomNeighbors.size() == 2 &&
atomRanks[begAtomNeighbors[0]] ==
atomRanks[begAtomNeighbors[1]]) ||
(endAtomNeighbors.size() == 2 &&
atomRanks[endAtomNeighbors[0]] ==
atomRanks[endAtomNeighbors[1]])) {
(*bondIt)->setStereo(Bond::STEREONONE);
(*bondIt)->getStereoAtoms().clear();
}
}
}
#endif
}
}
mol.setProp(common_properties::_StereochemDone, 1, true);
#if 0
std::cerr << "---\n";
mol.debugMol(std::cerr);
std::cerr << "<<<<<<<<<<<<<<<<\n";
#endif
}
// Find bonds than can be cis/trans in a molecule and mark them as
// Bond::STEREOANY.
void findPotentialStereoBonds(ROMol &mol, bool cleanIt) {
// FIX: The earlier thought was to provide an optional argument to ignore or
// consider
// double bonds in a ring. But I am removing this optional argument and
// ignoring ring bonds
// completely for now. This is because finding a potential stereo bond in a
// ring involves
// more than just checking the CIPranks for the neighbors - SP 05/04/04
// make this function callable multiple times
if ((mol.hasProp(common_properties::_BondsPotentialStereo)) && (!cleanIt)) {
return;
} else {
UINT_VECT ranks;
ranks.resize(mol.getNumAtoms());
bool cipDone = false;
ROMol::BondIterator bondIt;
for (bondIt = mol.beginBonds(); bondIt != mol.endBonds(); ++bondIt) {
if ((*bondIt)->getBondType() == Bond::DOUBLE &&
!(mol.getRingInfo()->numBondRings((*bondIt)->getIdx()))) {
// we are ignoring ring bonds here - read the FIX above
Bond *dblBond = *bondIt;
// We ignore bonds flagged as EITHERDOUBLE or STEREOANY which have
// stereo atoms set.
if (dblBond->getBondDir() == Bond::EITHERDOUBLE ||
(dblBond->getStereo() == Bond::STEREOANY &&
dblBond->getStereoAtoms().size() == 2)) {
continue;
}
// proceed only if we either want to clean the stereocode on this bond,
// if none is set on it yet, or it is STEREOANY and we need to find
// stereoatoms
if (cleanIt || dblBond->getStereo() == Bond::STEREONONE ||
dblBond->getStereo() == Bond::STEREOANY) {
dblBond->setStereo(Bond::STEREONONE);
const Atom *begAtom = dblBond->getBeginAtom(),
*endAtom = dblBond->getEndAtom();
// we're only going to handle 2 or three coordinate atoms:
if ((begAtom->getDegree() == 2 || begAtom->getDegree() == 3) &&
(endAtom->getDegree() == 2 || endAtom->getDegree() == 3)) {
// ------------------
// get the CIP ranking of each atom if we need it:
if (!cipDone) {
if (!begAtom->hasProp(common_properties::_CIPRank)) {
Chirality::assignAtomCIPRanks(mol, ranks);
} else {
// no need to recompute if we don't need to recompute. :-)
for (unsigned int ai = 0; ai < mol.getNumAtoms(); ++ai) {
ranks[ai] = mol.getAtomWithIdx(ai)->getProp<unsigned int>(
common_properties::_CIPRank);
}
}
cipDone = true;
}
// find the neighbors for the begin atom and the endAtom
UINT_VECT begAtomNeighbors, endAtomNeighbors;
bool checkDir = false;
bool includeAromatic = true;
Chirality::findAtomNeighborsHelper(mol, begAtom, dblBond,
begAtomNeighbors, checkDir,
includeAromatic);
Chirality::findAtomNeighborsHelper(mol, endAtom, dblBond,
endAtomNeighbors, checkDir,
includeAromatic);
if (begAtomNeighbors.size() > 0 && endAtomNeighbors.size() > 0) {
if ((begAtomNeighbors.size() == 2) &&
(endAtomNeighbors.size() == 2)) {
// if both of the atoms have 2 neighbors (other than the one
// connected
// by the double bond) and ....
#if 0
std::cerr << "Bond: " << dblBond->getIdx() << " "
<< begAtom->getIdx() << "=" << endAtom->getIdx()
<< std::endl;
std::cerr << " " << begAtomNeighbors[0] << "="
<< ranks[begAtomNeighbors[0]] << ":";
std::cerr << " " << begAtomNeighbors[1] << "="
<< ranks[begAtomNeighbors[1]] << std::endl;
std::cerr << " " << endAtomNeighbors[0] << "="
<< ranks[endAtomNeighbors[0]] << ":";
std::cerr << " " << endAtomNeighbors[1] << "="
<< ranks[endAtomNeighbors[1]] << std::endl;
#endif
if ((ranks[begAtomNeighbors[0]] !=
ranks[begAtomNeighbors[1]]) &&
(ranks[endAtomNeighbors[0]] !=
ranks[endAtomNeighbors[1]])) {
// the neighbors ranks are different at both the ends,
// this bond can be part of a cis/trans system
if (ranks[begAtomNeighbors[0]] > ranks[begAtomNeighbors[1]]) {
dblBond->getStereoAtoms().push_back(begAtomNeighbors[0]);
} else {
dblBond->getStereoAtoms().push_back(begAtomNeighbors[1]);
}
if (ranks[endAtomNeighbors[0]] > ranks[endAtomNeighbors[1]]) {
dblBond->getStereoAtoms().push_back(endAtomNeighbors[0]);
} else {
dblBond->getStereoAtoms().push_back(endAtomNeighbors[1]);
}
}
} else if (begAtomNeighbors.size() == 2) {
// if the begAtom has two neighbors and ....
if (ranks[begAtomNeighbors[0]] != ranks[begAtomNeighbors[1]]) {
// their ranks are different
if (ranks[begAtomNeighbors[0]] > ranks[begAtomNeighbors[1]]) {
dblBond->getStereoAtoms().push_back(begAtomNeighbors[0]);
} else {
dblBond->getStereoAtoms().push_back(begAtomNeighbors[1]);
}
dblBond->getStereoAtoms().push_back(endAtomNeighbors[0]);
}
} else if (endAtomNeighbors.size() == 2) {
// if the endAtom has two neighbors and ...
if (ranks[endAtomNeighbors[0]] != ranks[endAtomNeighbors[1]]) {
// their ranks are different
dblBond->getStereoAtoms().push_back(begAtomNeighbors[0]);
if (ranks[endAtomNeighbors[0]] > ranks[endAtomNeighbors[1]]) {
dblBond->getStereoAtoms().push_back(endAtomNeighbors[0]);
} else {
dblBond->getStereoAtoms().push_back(endAtomNeighbors[1]);
}
}
} else {
// end and beg atoms has only one neighbor each, it doesn't
// matter what the ranks are:
dblBond->getStereoAtoms().push_back(begAtomNeighbors[0]);
dblBond->getStereoAtoms().push_back(endAtomNeighbors[0]);
} // end of different number of neighbors on beg and end atoms
// mark this double bond as a potential stereo bond
if (!dblBond->getStereoAtoms().empty()) {
dblBond->setStereo(Bond::STEREOANY);
}
} // end of check that beg and end atoms have at least 1
// neighbor:
} // end of 2 and 3 coordinated atoms only
} // end of we want it or CIP code is not set
} // end of double bond
} // end of for loop over all bonds
mol.setProp(common_properties::_BondsPotentialStereo, 1, true);
}
}
// removes chirality markers from sp and sp2 hybridized centers:
void cleanupChirality(RWMol &mol) {
for (ROMol::AtomIterator atomIt = mol.beginAtoms(); atomIt != mol.endAtoms();
++atomIt) {
if ((*atomIt)->getChiralTag() != Atom::CHI_UNSPECIFIED &&
(*atomIt)->getHybridization() < Atom::SP3) {
(*atomIt)->setChiralTag(Atom::CHI_UNSPECIFIED);
}
}
}
void assignChiralTypesFrom3D(ROMol &mol, int confId, bool replaceExistingTags) {
const double ZERO_VOLUME_TOL = 0.1;
if (!mol.getNumConformers()) {
return;
}
const Conformer &conf = mol.getConformer(confId);
if (!conf.is3D()) {
return;
}
// if the molecule already has stereochemistry
// perceived, remove the flags that indicate
// this... what we're about to do will require
// that we go again.
if (mol.hasProp(common_properties::_StereochemDone)) {
mol.clearProp(common_properties::_StereochemDone);
}
for (ROMol::AtomIterator atomIt = mol.beginAtoms(); atomIt != mol.endAtoms();
++atomIt) {
Atom *atom = *atomIt;
// if we aren't replacing existing tags and the atom is already tagged,
// punt:
if (!replaceExistingTags && atom->getChiralTag() != Atom::CHI_UNSPECIFIED) {
continue;
}
atom->setChiralTag(Atom::CHI_UNSPECIFIED);
// additional reasons to skip the atom:
if (atom->getDegree() < 3 || atom->getTotalDegree() > 4) {
// not enough explicit neighbors or too many total neighbors
continue;
} else {
int anum = atom->getAtomicNum();
if (anum != 16 && anum != 34 && // S or Se are special
// (just using the InChI list for now)
(atom->getTotalDegree() != 4 || // not enough total neighbors
atom->getTotalNumHs(true) > 1)) {
continue;
}
}
const RDGeom::Point3D &p0 = conf.getAtomPos(atom->getIdx());
ROMol::ADJ_ITER nbrIdx, endNbrs;
boost::tie(nbrIdx, endNbrs) = mol.getAtomNeighbors(atom);
const RDGeom::Point3D &p1 = conf.getAtomPos(*nbrIdx);
++nbrIdx;
const RDGeom::Point3D &p2 = conf.getAtomPos(*nbrIdx);
++nbrIdx;
const RDGeom::Point3D &p3 = conf.getAtomPos(*nbrIdx);
RDGeom::Point3D v1 = p1 - p0;
RDGeom::Point3D v2 = p2 - p0;
RDGeom::Point3D v3 = p3 - p0;
double chiralVol = v1.dotProduct(v2.crossProduct(v3));
if (chiralVol < -ZERO_VOLUME_TOL) {
atom->setChiralTag(Atom::CHI_TETRAHEDRAL_CW);
} else if (chiralVol > ZERO_VOLUME_TOL) {
atom->setChiralTag(Atom::CHI_TETRAHEDRAL_CCW);
} else {
atom->setChiralTag(Atom::CHI_UNSPECIFIED);
}
}
}
void assignChiralTypesFromMolParity(ROMol &mol, bool replaceExistingTags) {
static const std::vector<Atom::ChiralType> chiralTypeVect{
Atom::CHI_TETRAHEDRAL_CW, Atom::CHI_TETRAHEDRAL_CCW};
// if the molecule already has stereochemistry
// perceived, remove the flags that indicate
// this... what we're about to do will require
// that we go again.
if (mol.hasProp(common_properties::_StereochemDone)) {
mol.clearProp(common_properties::_StereochemDone);
}
// Atom-based parity
// Number the atoms surrounding the stereo center with 1, 2, 3, and 4
// in order of increasing atom number (position in the atom block)
// (an implicit hydrogen should be considered the highest numbered atom).
// View the center from a position such that the bond connecting the
// highest-numbered atom (4) projects behind the plane formed by
// atoms 1, 2, and 3.
//
// Parity 1 (CW) Parity 2 (CCW)
// 3 1 3 2
// \ / \ /
// | |
// 2 1
//
for (auto atom : mol.atoms()) {
// if we aren't replacing existing tags and the atom is already tagged,
// punt:
if (!replaceExistingTags && atom->getChiralTag() != Atom::CHI_UNSPECIFIED) {
continue;
}
int parity = 0;
atom->getPropIfPresent(common_properties::molParity, parity);
if (parity <= 0 || parity > 2 || atom->getDegree() < 3) {
atom->setChiralTag(Atom::CHI_UNSPECIFIED);
continue;
}
// if we are here, parity was 1 (CW) or 2 (CCW)
// now we set parity 0 to be CW and 1 to be CCW
--parity;
RDKit::ROMol::OBOND_ITER_PAIR nbrBonds = mol.getAtomBonds(atom);
INT_LIST nbrBondIdxList;
std::transform(
nbrBonds.first, nbrBonds.second, std::back_inserter(nbrBondIdxList),
[mol](const ROMol::edge_descriptor &e) { return mol[e]->getIdx(); });
unsigned int atomIdx = atom->getIdx();
nbrBondIdxList.sort([mol, atomIdx](const int ai, const int bi) {
return (mol.getBondWithIdx(ai)->getOtherAtomIdx(atomIdx) <
mol.getBondWithIdx(bi)->getOtherAtomIdx(atomIdx));
});
int nSwaps = atom->getPerturbationOrder(nbrBondIdxList);
if (nSwaps % 2) {
parity = 1 - parity;
}
atom->setChiralTag(chiralTypeVect[parity]);
if (atom->getImplicitValence() == -1) {
atom->calcExplicitValence(false);
atom->calcImplicitValence(false);
}
// within the RD representation, if a three-coordinate atom
// is chiral and has an implicit H, that H needs to be made explicit:
if (atom->getDegree() == 3 && !atom->getNumExplicitHs() &&
atom->getNumImplicitHs() == 1) {
atom->setNumExplicitHs(1);
// recalculated number of implicit Hs:
atom->updatePropertyCache();
}
}
}
void setDoubleBondNeighborDirections(ROMol &mol, const Conformer *conf) {
// used to store the number of single bonds a given
// single bond is adjacent to
std::vector<unsigned int> singleBondCounts(mol.getNumBonds(), 0);
std::vector<Bond *> bondsInPlay;
// keeps track of which single bonds are adjacent to each double bond:
VECT_INT_VECT dblBondNbrs(mol.getNumBonds());
// keeps track of which double bonds are adjacent to each single bond:
VECT_INT_VECT singleBondNbrs(mol.getNumBonds());
// keeps track of which single bonds need a dir set and which double bonds
// need to have their neighbors' dirs set
boost::dynamic_bitset<> needsDir(mol.getNumBonds());
// find double bonds that should be considered for
// stereochemistry
// NOTE that we are explicitly excluding double bonds in rings
// with this test.
bool resetRings = false;
if (!mol.getRingInfo()->isInitialized()) {
resetRings = true;
MolOps::fastFindRings(mol);
}
for (RWMol::BondIterator bondIt = mol.beginBonds(); bondIt != mol.endBonds();
++bondIt) {
if (isBondCandidateForStereo(*bondIt)) {
const Atom *a1 = (*bondIt)->getBeginAtom();
const Atom *a2 = (*bondIt)->getEndAtom();
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(a1);
while (beg != end) {
const Bond *nbrBond = mol[*beg];
if (nbrBond->getBondType() == Bond::SINGLE ||
nbrBond->getBondType() == Bond::AROMATIC) {
singleBondCounts[nbrBond->getIdx()] += 1;
auto nbrDir = nbrBond->getBondDir();
if (nbrDir == Bond::BondDir::NONE ||
nbrDir == Bond::BondDir::ENDDOWNRIGHT ||
nbrDir == Bond::BondDir::ENDUPRIGHT) {
needsDir[nbrBond->getIdx()] = 1;
}
needsDir[(*bondIt)->getIdx()] = 1;
dblBondNbrs[(*bondIt)->getIdx()].push_back(nbrBond->getIdx());
// the search may seem inefficient, but these vectors are going to
// be at most 2 long (with very few exceptions). It's just not worth
// using a different data structure
if (std::find(singleBondNbrs[nbrBond->getIdx()].begin(),
singleBondNbrs[nbrBond->getIdx()].end(),
(*bondIt)->getIdx()) ==
singleBondNbrs[nbrBond->getIdx()].end()) {
singleBondNbrs[nbrBond->getIdx()].push_back((*bondIt)->getIdx());
}
}
++beg;
}
boost::tie(beg, end) = mol.getAtomBonds(a2);
while (beg != end) {
const Bond *nbrBond = mol[*beg];
if (nbrBond->getBondType() == Bond::SINGLE ||
nbrBond->getBondType() == Bond::AROMATIC) {
singleBondCounts[nbrBond->getIdx()] += 1;
auto nbrDir = nbrBond->getBondDir();
if (nbrDir == Bond::BondDir::NONE ||
nbrDir == Bond::BondDir::ENDDOWNRIGHT ||
nbrDir == Bond::BondDir::ENDUPRIGHT) {
needsDir[nbrBond->getIdx()] = 1;
}
needsDir[(*bondIt)->getIdx()] = 1;
dblBondNbrs[(*bondIt)->getIdx()].push_back(nbrBond->getIdx());
// the search may seem inefficient, but these vectors are going to
// be at most 2 long (with very few exceptions). It's just not worth
// using a different data structure
if (std::find(singleBondNbrs[nbrBond->getIdx()].begin(),
singleBondNbrs[nbrBond->getIdx()].end(),
(*bondIt)->getIdx()) ==
singleBondNbrs[nbrBond->getIdx()].end()) {
singleBondNbrs[nbrBond->getIdx()].push_back((*bondIt)->getIdx());
}
}
++beg;
}
bondsInPlay.push_back(*bondIt);
}
}
if (!bondsInPlay.size()) {
if (resetRings) {
mol.getRingInfo()->reset();
}
return;
}
// order the double bonds based on the singleBondCounts of their neighbors:
std::vector<std::pair<unsigned int, Bond *>> orderedBondsInPlay;
for (auto dblBond : bondsInPlay) {
unsigned int countHere =
std::accumulate(dblBondNbrs[dblBond->getIdx()].begin(),
dblBondNbrs[dblBond->getIdx()].end(), 0);
// and favor double bonds that are *not* in rings. The combination of
// using the sum above (instead of the max) and this ring-membershipt test
// seem to fix sf.net issue 3009836
if (!(mol.getRingInfo()->numBondRings(dblBond->getIdx()))) {
countHere *= 10;
}
orderedBondsInPlay.push_back(std::make_pair(countHere, dblBond));
}
std::sort(orderedBondsInPlay.begin(), orderedBondsInPlay.end());
// oof, now loop over the double bonds in that order and
// update their neighbor directionalities:
std::vector<std::pair<unsigned int, Bond *>>::reverse_iterator pairIter;
for (pairIter = orderedBondsInPlay.rbegin();
pairIter != orderedBondsInPlay.rend(); ++pairIter) {
// std::cerr << "RESET?: " << pairIter->second->getIdx() << " "
// << pairIter->second->getStereo() << std::endl;
updateDoubleBondNeighbors(mol, pairIter->second, conf, needsDir,
singleBondCounts, singleBondNbrs);
}
if (resetRings) {
mol.getRingInfo()->reset();
}
}
void detectBondStereochemistry(ROMol &mol, int confId) {
if (!mol.getNumConformers()) {
return;
}
const Conformer &conf = mol.getConformer(confId);
setDoubleBondNeighborDirections(mol, &conf);
}
void setBondStereoFromDirections(ROMol &mol) {
for (Bond *bond : mol.bonds()) {
if (bond->getBondType() == Bond::DOUBLE) {
const Atom *stereoBondBeginAtom = bond->getBeginAtom();
const Atom *stereoBondEndAtom = bond->getEndAtom();
const Bond *directedBondAtBegin =
Chirality::getNeighboringDirectedBond(mol, stereoBondBeginAtom);
const Bond *directedBondAtEnd =
Chirality::getNeighboringDirectedBond(mol, stereoBondEndAtom);
if (directedBondAtBegin != nullptr && directedBondAtEnd != nullptr) {
unsigned beginSideStereoAtom =
directedBondAtBegin->getOtherAtomIdx(stereoBondBeginAtom->getIdx());
unsigned endSideStereoAtom =
directedBondAtEnd->getOtherAtomIdx(stereoBondEndAtom->getIdx());
bond->setStereoAtoms(beginSideStereoAtom, endSideStereoAtom);
auto beginSideBondDirection = directedBondAtBegin->getBondDir();
if (directedBondAtBegin->getBeginAtom() == stereoBondBeginAtom) {
beginSideBondDirection = getOppositeBondDir(beginSideBondDirection);
}
auto endSideBondDirection = directedBondAtEnd->getBondDir();
if (directedBondAtEnd->getEndAtom() == stereoBondEndAtom) {
endSideBondDirection = getOppositeBondDir(endSideBondDirection);
}
if (beginSideBondDirection == endSideBondDirection) {
bond->setStereo(Bond::STEREOTRANS);
} else {
bond->setStereo(Bond::STEREOCIS);
}
}
}
}
}
void assignStereochemistryFrom3D(ROMol &mol, int confId,
bool replaceExistingTags) {
if (!mol.getNumConformers() || !mol.getConformer(confId).is3D()) {
return;
}
detectBondStereochemistry(mol, confId);
assignChiralTypesFrom3D(mol, confId, replaceExistingTags);
bool force = true;
bool flagPossibleStereoCenters = true;
assignStereochemistry(mol, replaceExistingTags, force,
flagPossibleStereoCenters);
}
void assignChiralTypesFromBondDirs(ROMol &mol, const int confId,
const bool replaceExistingTags) {
if (!mol.getNumConformers()) {
return;
}
auto conf = mol.getConformer(confId);
boost::dynamic_bitset<> atomsSet(mol.getNumAtoms(), 0);
for (auto &bond : mol.bonds()) {
const Bond::BondDir dir = bond->getBondDir();
if (dir != Bond::UNKNOWN) {
// the bond is marked as chiral:
if (dir == Bond::BEGINWEDGE || dir == Bond::BEGINDASH) {
Atom *atom = bond->getBeginAtom();
if (atomsSet[atom->getIdx()] ||
(!replaceExistingTags &&
atom->getChiralTag() != Atom::CHI_UNSPECIFIED)) {
continue;
}
if (atom->getImplicitValence() == -1) {
atom->calcExplicitValence(false);
atom->calcImplicitValence(false);
}
Atom::ChiralType code = atomChiralTypeFromBondDir(mol, bond, &conf);
if (code != Atom::ChiralType::CHI_UNSPECIFIED) {
atomsSet.set(atom->getIdx());
// std::cerr << "atom " << atom->getIdx() << " code " << code
// << " from bond " << bond->getIdx() << std::endl;
}
atom->setChiralTag(code);
// within the RD representation, if a three-coordinate atom
// is chiral and has an implicit H, that H needs to be made explicit:
if (atom->getDegree() == 3 && !atom->getNumExplicitHs() &&
atom->getNumImplicitHs() == 1) {
atom->setNumExplicitHs(1);
// recalculated number of implicit Hs:
atom->updatePropertyCache();
}
}
}
}
}
void removeStereochemistry(ROMol &mol) {
if (mol.hasProp(common_properties::_StereochemDone)) {
mol.clearProp(common_properties::_StereochemDone);
}
for (ROMol::AtomIterator atIt = mol.beginAtoms(); atIt != mol.endAtoms();
++atIt) {
(*atIt)->setChiralTag(Atom::CHI_UNSPECIFIED);
if ((*atIt)->hasProp(common_properties::_CIPCode)) {
(*atIt)->clearProp(common_properties::_CIPCode);
}
if ((*atIt)->hasProp(common_properties::_CIPRank)) {
(*atIt)->clearProp(common_properties::_CIPRank);
}
}
for (ROMol::BondIterator bondIt = mol.beginBonds(); bondIt != mol.endBonds();
++bondIt) {
if ((*bondIt)->getBondType() == Bond::DOUBLE) {
(*bondIt)->setStereo(Bond::STEREONONE);
(*bondIt)->getStereoAtoms().clear();
} else if ((*bondIt)->getBondType() == Bond::SINGLE) {
(*bondIt)->setBondDir(Bond::NONE);
}
}
std::vector<StereoGroup> sgs;
static_cast<RWMol &>(mol).setStereoGroups(std::move(sgs));
}
} // end of namespace MolOps
} // end of namespace RDKit
Reference in New Issue View Git Blame Copy Permalink