pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-03 21:44:30 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
Release_2021_03
rdkit/Code/GraphMol/Canon.cpp
jfkonecn 435a6b64b5 fixes #3967 (#4190) 
* fixes #3967

* modernize test

Co-authored-by: greg landrum <greg.landrum@gmail.com>
2021-06-09 06:12:31 +02:00

1259 lines
49 KiB
C++
Raw Permalink Blame History

//
// Copyright (C) 2001-2020 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 <GraphMol/Canon.h>
#include <GraphMol/SmilesParse/SmilesParseOps.h>
#include <GraphMol/RDKitQueries.h>
#include <RDGeneral/Exceptions.h>
#include <RDGeneral/hash/hash.hpp>
#include <RDGeneral/utils.h>
#include <algorithm>
namespace RDKit {
namespace Canon {
namespace {
bool isUnsaturated(const Atom *atom, const ROMol &mol) {
for (const auto &bndItr :
boost::make_iterator_range(mol.getAtomBonds(atom))) {
// can't just check for single bonds, because dative bonds also have an
// order of 1
if (mol[bndItr]->getBondTypeAsDouble() > 1) {
return true;
}
}
return false;
}
bool hasSingleHQuery(const Atom::QUERYATOM_QUERY *q) {
// list queries are series of nested ors of AtomAtomicNum queries
PRECONDITION(q, "bad query");
bool res = false;
std::string descr = q->getDescription();
if (descr == "AtomAnd") {
for (auto cIt = q->beginChildren(); cIt != q->endChildren(); ++cIt) {
std::string descr = (*cIt)->getDescription();
if (descr == "AtomHCount") {
if (!(*cIt)->getNegation() &&
((ATOM_EQUALS_QUERY *)(*cIt).get())->getVal() == 1) {
return true;
}
return false;
} else if (descr == "AtomAnd") {
res = hasSingleHQuery((*cIt).get());
if (res) {
return true;
}
}
}
}
return res;
}
bool atomHasFourthValence(const Atom *atom) {
if (atom->getNumExplicitHs() == 1 || atom->getImplicitValence() == 1) {
return true;
}
if (atom->hasQuery()) {
// the SMARTS [C@@H] produces an atom with a H query, but we also
// need to treat this like an explicit H for chirality purposes
// This was Github #1489
return hasSingleHQuery(atom->getQuery());
}
return false;
}
} // end of anonymous namespace
bool chiralAtomNeedsTagInversion(const RDKit::ROMol &mol,
const RDKit::Atom *atom, bool isAtomFirst,
size_t numClosures) {
PRECONDITION(atom, "bad atom");
return atom->getDegree() == 3 &&
((isAtomFirst && atom->getNumExplicitHs() == 1) ||
(!atomHasFourthValence(atom) && numClosures == 1 &&
!isUnsaturated(atom, mol)));
}
struct _possibleCompare
: public std::binary_function<PossibleType, PossibleType, bool> {
bool operator()(const PossibleType &arg1, const PossibleType &arg2) const {
return (arg1.get<0>() < arg2.get<0>());
}
};
bool checkBondsInSameBranch(MolStack &molStack, Bond *dblBnd, Bond *dirBnd) {
bool seenDblBond = false;
int branchCounter = 0;
for (const auto &item : molStack) {
switch (item.type) {
case MOL_STACK_BOND:
if (item.obj.bond == dirBnd || item.obj.bond == dblBnd) {
if (seenDblBond) {
return branchCounter == 0;
} else {
seenDblBond = true;
}
}
break;
case MOL_STACK_BRANCH_OPEN:
if (seenDblBond) {
++branchCounter;
}
break;
case MOL_STACK_BRANCH_CLOSE:
if (seenDblBond) {
--branchCounter;
}
break;
default:
break;
}
}
// We should not ever hit this. But if we do, returning false
// causes the same behavior as before this patch.
return false;
}
void switchBondDir(Bond *bond) {
PRECONDITION(bond, "bad bond");
PRECONDITION(bond->getBondType() == Bond::SINGLE || bond->getIsAromatic(),
"bad bond type");
switch (bond->getBondDir()) {
case Bond::ENDUPRIGHT:
bond->setBondDir(Bond::ENDDOWNRIGHT);
break;
case Bond::ENDDOWNRIGHT:
bond->setBondDir(Bond::ENDUPRIGHT);
break;
default:
break;
}
}
// FIX: this may only be of interest from the SmilesWriter, should we
// move it there?
//
//
void canonicalizeDoubleBond(Bond *dblBond, UINT_VECT &bondVisitOrders,
UINT_VECT &atomVisitOrders,
UINT_VECT &bondDirCounts, UINT_VECT &atomDirCounts,
MolStack &molStack) {
PRECONDITION(dblBond, "bad bond");
PRECONDITION(dblBond->getBondType() == Bond::DOUBLE, "bad bond order");
PRECONDITION(dblBond->getStereo() > Bond::STEREOANY, "bad bond stereo");
PRECONDITION(dblBond->getStereoAtoms().size() >= 2, "bad bond stereo atoms");
PRECONDITION(atomVisitOrders[dblBond->getBeginAtomIdx()] > 0 ||
atomVisitOrders[dblBond->getEndAtomIdx()] > 0,
"neither end atom traversed");
Atom *atom1 = dblBond->getBeginAtom();
Atom *atom2 = dblBond->getEndAtom();
// we only worry about double bonds that begin and end at atoms
// of degree 2 or 3:
if ((atom1->getDegree() != 2 && atom1->getDegree() != 3) ||
(atom2->getDegree() != 2 && atom2->getDegree() != 3)) {
return;
}
// ensure that atom1 is the lower numbered atom of the double bond (the one
// traversed first)
if (atomVisitOrders[dblBond->getBeginAtomIdx()] >=
atomVisitOrders[dblBond->getEndAtomIdx()]) {
std::swap(atom1, atom2);
}
Bond *firstFromAtom1 = nullptr, *secondFromAtom1 = nullptr;
Bond *firstFromAtom2 = nullptr, *secondFromAtom2 = nullptr;
ROMol &mol = dblBond->getOwningMol();
auto firstVisitOrder = mol.getNumBonds() + 1;
ROMol::OBOND_ITER_PAIR atomBonds;
// -------------------------------------------------------
// find the lowest visit order bonds from each end and determine
// if anything is already constraining our choice of directions:
bool dir1Set = false, dir2Set = false;
for (const auto &bndItr :
boost::make_iterator_range(mol.getAtomBonds(atom1))) {
auto bond = mol[bndItr];
if (bond != dblBond) {
auto bondIdx = bond->getIdx();
if (bondDirCounts[bondIdx] > 0) {
dir1Set = true;
}
if (!firstFromAtom1 || bondVisitOrders[bondIdx] < firstVisitOrder) {
if (firstFromAtom1) {
secondFromAtom1 = firstFromAtom1;
}
firstFromAtom1 = bond;
firstVisitOrder = bondVisitOrders[bondIdx];
} else {
secondFromAtom1 = bond;
}
}
}
firstVisitOrder = mol.getNumBonds() + 1;
for (const auto &bndItr :
boost::make_iterator_range(mol.getAtomBonds(atom2))) {
auto bond = mol[bndItr];
if (bond != dblBond) {
auto bondIdx = bond->getIdx();
if (bondDirCounts[bondIdx] > 0) {
dir2Set = true;
}
if (!firstFromAtom2 || bondVisitOrders[bondIdx] < firstVisitOrder) {
if (firstFromAtom2) {
secondFromAtom2 = firstFromAtom2;
}
firstFromAtom2 = bond;
firstVisitOrder = bondVisitOrders[bondIdx];
} else {
secondFromAtom2 = bond;
}
}
}
// make sure we found everything we need to find:
CHECK_INVARIANT(firstFromAtom1, "could not find atom1");
CHECK_INVARIANT(firstFromAtom2, "could not find atom2");
CHECK_INVARIANT(atom1->getDegree() == 2 || secondFromAtom1,
"inconsistency at atom1");
CHECK_INVARIANT(atom2->getDegree() == 2 || secondFromAtom2,
"inconsistency at atom2");
bool setFromBond1 = true;
Bond::BondDir atom1Dir = Bond::NONE;
Bond::BondDir atom2Dir = Bond::NONE;
Bond *atom1ControllingBond = firstFromAtom1;
Bond *atom2ControllingBond = firstFromAtom2;
if (!dir1Set && !dir2Set) {
// ----------------------------------
// nothing has touched our bonds so far, so set the
// directions to "arbitrary" values:
// the bond we came in on becomes ENDUPRIGHT:
atom1Dir = Bond::ENDUPRIGHT;
firstFromAtom1->setBondDir(atom1Dir);
bondDirCounts[firstFromAtom1->getIdx()] += 1;
atomDirCounts[atom1->getIdx()] += 1;
} else if (!dir2Set) {
// at least one of the bonds on atom1 has its directionality set already:
if (bondDirCounts[firstFromAtom1->getIdx()] > 0) {
// The first bond's direction has been set at some earlier point:
atom1Dir = firstFromAtom1->getBondDir();
bondDirCounts[firstFromAtom1->getIdx()] += 1;
atomDirCounts[atom1->getIdx()] += 1;
if (secondFromAtom1) {
// both bonds have their directionalities set, make sure
// they are compatible:
if (firstFromAtom1->getBondDir() == secondFromAtom1->getBondDir() &&
bondDirCounts[firstFromAtom2->getIdx()]) {
CHECK_INVARIANT(
((firstFromAtom1->getBeginAtomIdx() == atom1->getIdx()) ^
(secondFromAtom1->getBeginAtomIdx() == atom1->getIdx())),
"inconsistent state");
}
}
} else {
// the second bond must be present and setting the direction:
CHECK_INVARIANT(secondFromAtom1, "inconsistent state");
CHECK_INVARIANT(bondDirCounts[secondFromAtom1->getIdx()] > 0,
"inconsistent state");
// It must be the second bond setting the direction.
// This happens when the bond dir is set in a branch:
// v- this double bond
// CC(/C=P/N)=N/O
// ^- the second bond sets the direction
// or when the first bond is a ring closure from an
// earlier traversed atom:
// v- this double bond
// NC1=NOC/C1=N\O
// ^- this closure ends up being the first bond,
// and it does not set the direction.
//
// This addresses parts of Issue 185 and sf.net Issue 1842174
//
atom1Dir = secondFromAtom1->getBondDir();
firstFromAtom1->setBondDir(atom1Dir);
bondDirCounts[firstFromAtom1->getIdx()] += 1;
atomDirCounts[atom1->getIdx()] += 2;
atom1ControllingBond = secondFromAtom1;
}
} else {
// dir2 has been set, and dir1 hasn't: we're dealing with a stereochem
// specification on a ring double bond:
setFromBond1 = false;
// at least one of the bonds on atom2 has its directionality set already:
if (bondDirCounts[firstFromAtom2->getIdx()] > 0) {
// The second bond's direction has been set at some earlier point:
atom2Dir = firstFromAtom2->getBondDir();
bondDirCounts[firstFromAtom2->getIdx()] += 1;
atomDirCounts[atom2->getIdx()] += 1;
if (secondFromAtom2) {
// both bonds have their directionalities set, make sure
// they are compatible:
if (firstFromAtom2->getBondDir() == secondFromAtom2->getBondDir() &&
bondDirCounts[firstFromAtom1->getIdx()]) {
CHECK_INVARIANT(
((firstFromAtom2->getBeginAtomIdx() == atom2->getIdx()) ^
(secondFromAtom2->getBeginAtomIdx() == atom2->getIdx())),
"inconsistent state");
}
}
} else {
// the second bond must be present and setting the direction:
CHECK_INVARIANT(secondFromAtom2, "inconsistent state");
CHECK_INVARIANT(bondDirCounts[secondFromAtom2->getIdx()] > 0,
"inconsistent state");
// It must be the second bond setting the direction.
// This happens when the bond dir is set in a branch:
// v- this double bond
// CC(/C=P/N)=N/O
// ^- the second bond sets the direction
// or when the first bond is a ring closure from an
// earlier traversed atom:
// v- this double bond
// NC1=NOC/C1=N\O
// ^- this closure ends up being the first bond,
// and it does not set the direction.
//
// This addresses parts of Issue 185 and sf.net Issue 1842174
//
atom2Dir = secondFromAtom2->getBondDir();
firstFromAtom2->setBondDir(atom2Dir);
bondDirCounts[firstFromAtom2->getIdx()] += 1;
atomDirCounts[atom2->getIdx()] += 2;
atom2ControllingBond = secondFromAtom2;
}
// CHECK_INVARIANT(0,"ring stereochemistry not handled");
} // end of the ring stereochemistry if
// now set the directionality on the other side:
if (setFromBond1) {
if (dblBond->getStereo() == Bond::STEREOE ||
dblBond->getStereo() == Bond::STEREOTRANS) {
atom2Dir = atom1Dir;
} else if (dblBond->getStereo() == Bond::STEREOZ ||
dblBond->getStereo() == Bond::STEREOCIS) {
atom2Dir = (atom1Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
}
CHECK_INVARIANT(atom2Dir != Bond::NONE, "stereo not set");
// If we're not looking at the bonds used to determine the
// stereochemistry, we need to flip the setting on the other bond:
const INT_VECT &stereoAtoms = dblBond->getStereoAtoms();
auto isClosingRingBond = [](Bond *bond) {
if (bond == nullptr) {
return false;
}
auto beginIdx = bond->getBeginAtomIdx();
auto endIdx = bond->getEndAtomIdx();
return beginIdx > endIdx &&
beginIdx - endIdx > 1 &&
bond->hasProp(common_properties::_TraversalRingClosureBond);
};
auto isFlipped = false;
if (atom1->getDegree() == 3 &&
std::find(stereoAtoms.begin(), stereoAtoms.end(),
static_cast<int>(atom1ControllingBond->getOtherAtomIdx(
atom1->getIdx()))) == stereoAtoms.end()) {
isFlipped = true;
atom2Dir = (atom2Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
}
// std::cerr<<" 0 set bond 2: "<<firstFromAtom2->getIdx()<<"
// "<<atom2Dir<<std::endl;
if (atom2->getDegree() == 3 &&
std::find(stereoAtoms.begin(), stereoAtoms.end(),
static_cast<int>(firstFromAtom2->getOtherAtomIdx(
atom2->getIdx()))) == stereoAtoms.end()) {
isFlipped = true;
atom2Dir = (atom2Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
}
if (!isFlipped &&
(isClosingRingBond(dblBond) ||
(isClosingRingBond(secondFromAtom1) &&
!secondFromAtom1->getIsAromatic() &&
secondFromAtom1->getBondDir() != Bond::NONE))
) {
atom2Dir = (atom2Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
}
// std::cerr<<" 1 set bond 2: "<<firstFromAtom2->getIdx()<<"
// "<<atom2Dir<<std::endl;
firstFromAtom2->setBondDir(atom2Dir);
bondDirCounts[firstFromAtom2->getIdx()] += 1;
atomDirCounts[atom2->getIdx()] += 1;
} else {
// we come before a ring closure:
if (dblBond->getStereo() == Bond::STEREOZ ||
dblBond->getStereo() == Bond::STEREOCIS) {
atom1Dir = atom2Dir;
} else if (dblBond->getStereo() == Bond::STEREOE ||
dblBond->getStereo() == Bond::STEREOTRANS) {
atom1Dir = (atom2Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
}
CHECK_INVARIANT(atom1Dir != Bond::NONE, "stereo not set");
// If we're not looking at the bonds used to determine the
// stereochemistry, we need to flip the setting on the other bond:
const INT_VECT &stereoAtoms = dblBond->getStereoAtoms();
if (atom2->getDegree() == 3 &&
std::find(stereoAtoms.begin(), stereoAtoms.end(),
static_cast<int>(atom2ControllingBond->getOtherAtomIdx(
atom2->getIdx()))) == stereoAtoms.end()) {
// std::cerr<<"flip 1"<<std::endl;
atom1Dir = (atom1Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
}
if (atom1->getDegree() == 3 &&
std::find(stereoAtoms.begin(), stereoAtoms.end(),
static_cast<int>(firstFromAtom1->getOtherAtomIdx(
atom1->getIdx()))) == stereoAtoms.end()) {
// std::cerr<<"flip 2"<<std::endl;
atom1Dir = (atom1Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
}
firstFromAtom1->setBondDir(atom1Dir);
switchBondDir(firstFromAtom1);
bondDirCounts[firstFromAtom1->getIdx()] += 1;
atomDirCounts[atom1->getIdx()] += 1;
}
// -----------------------------------
//
// Check if there are other bonds from atoms 1 and 2 that need
// to have their directionalities set:
///
if (atom1->getDegree() == 3) {
if (!bondDirCounts[secondFromAtom1->getIdx()]) {
// This bond (the second bond from the starting atom of the double bond)
// is a special case. It's going to appear in a branch in the smiles:
// X\C(\Y)=C/Z
// ^
// |- here
// so it actually needs to go down with the *same* direction as the
// bond that's already been set (because "pulling the bond out of the
// branch" reverses its direction).
// A quick example. This SMILES:
// F/C(\Cl)=C/F
// is *wrong*. This is the correct form:
// F/C(/Cl)=C/F
// So, since we want this bond to have the opposite direction to the
// other one, we put it in with the same direction.
// This was Issue 183
// UNLESS the bond is not in a branch (in the smiles) (e.g. firstFromAtom1
// branches off a cycle, and secondFromAtom1 shows up at the end of the
// cycle). This was Github Issue #2023, see it for an example.
if (checkBondsInSameBranch(molStack, dblBond, secondFromAtom1)) {
auto otherDir = (firstFromAtom1->getBondDir() == Bond::ENDUPRIGHT)
? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
secondFromAtom1->setBondDir(otherDir);
} else {
secondFromAtom1->setBondDir(firstFromAtom1->getBondDir());
}
}
bondDirCounts[secondFromAtom1->getIdx()] += 1;
atomDirCounts[atom1->getIdx()] += 1;
}
if (atom2->getDegree() == 3) {
if (!bondDirCounts[secondFromAtom2->getIdx()]) {
// Here we set the bond direction to be opposite the other one (since
// both come after the atom connected to the double bond).
Bond::BondDir otherDir;
if (!secondFromAtom2->hasProp(
common_properties::_TraversalRingClosureBond)) {
otherDir = (firstFromAtom2->getBondDir() == Bond::ENDUPRIGHT)
? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
} else {
// another one those irritating little reversal things due to
// ring closures
otherDir = firstFromAtom2->getBondDir();
}
secondFromAtom2->setBondDir(otherDir);
}
bondDirCounts[secondFromAtom2->getIdx()] += 1;
atomDirCounts[atom2->getIdx()] += 1;
// std::cerr<<" other: "<<secondFromAtom2->getIdx()<<"
// "<<otherDir<<std::endl;
}
if (setFromBond1) {
// This is an odd case... The bonds off the beginning atom are
// after the start atom in the traversal stack. These need to
// have their directions reversed. An example SMILES (unlikely
// to actually traverse this way is:
// C(=C/O)/F or C(/F)=C/O
// That bond is Z, without the reversal, this would come out:
// C(=C/O)\F or C(\F)=C/O
// which is E.
//
// In the case of three-coordinate atoms, we don't need to flip
// the second bond because the Issue 183 fix (above) already got
// that one.
//
// This was Issue 191 and continued into sf.net issue 1842174
if (bondVisitOrders[atom1ControllingBond->getIdx()] >
atomVisitOrders[atom1->getIdx()]) {
if (bondDirCounts[atom1ControllingBond->getIdx()] == 1) {
if (!atom1ControllingBond->hasProp(
common_properties::_TraversalRingClosureBond)) {
// std::cerr<<" switcheroo 1"<<std::endl;
switchBondDir(atom1ControllingBond);
}
} else if (bondDirCounts[firstFromAtom2->getIdx()] == 1) {
// the controlling bond at atom1 is being set by someone else, flip the
// direction
// on the atom2 bond instead:
// std::cerr<<" switcheroo 2"<<std::endl;
switchBondDir(firstFromAtom2);
if (secondFromAtom2 && bondDirCounts[secondFromAtom2->getIdx()] >= 1) {
switchBondDir(secondFromAtom2);
}
}
}
}
// something to watch out for here. For this molecule and traversal order:
// 0 1 2 3 4 5 6 7 8 <- atom numbers
// C/C=C/C(/N=C/C)=C/C
// ^ ^
// |--|-- these two bonds must match in direction or the SMILES
// is inconsistent (according to Daylight, Marvin does ok with
// it)
// That means that the direction of the bond from atom 3->4 needs to be set
// when the bond from 2->3 is set.
// Issue2023: But only if 3->4 doesn't have a direction yet?
//
// I believe we only need to worry about this for the bonds from atom2.
const Atom *atom3 = firstFromAtom2->getOtherAtom(atom2);
if (atom3->getDegree() == 3) {
Bond *otherAtom3Bond = nullptr;
bool dblBondPresent = false;
atomBonds = mol.getAtomBonds(atom3);
while (atomBonds.first != atomBonds.second) {
Bond *tbond = mol[*atomBonds.first];
if (tbond->getBondType() == Bond::DOUBLE &&
tbond->getStereo() > Bond::STEREOANY) {
dblBondPresent = true;
} else if ((tbond->getBondType() == Bond::SINGLE) &&
(tbond != firstFromAtom2)) {
otherAtom3Bond = tbond;
}
atomBonds.first++;
}
if (dblBondPresent && otherAtom3Bond &&
otherAtom3Bond->getBondDir() == Bond::NONE) {
// std::cerr<<"set!"<<std::endl;
otherAtom3Bond->setBondDir(firstFromAtom2->getBondDir());
bondDirCounts[otherAtom3Bond->getIdx()] += 1;
atomDirCounts[atom3->getIdx()] += 1;
}
}
}
// finds cycles
void dfsFindCycles(ROMol &mol, int atomIdx, int inBondIdx,
std::vector<AtomColors> &colors, const UINT_VECT &ranks,
UINT_VECT &atomOrders, VECT_INT_VECT &atomRingClosures,
const boost::dynamic_bitset<> *bondsInPlay,
const std::vector<std::string> *bondSymbols, bool doRandom) {
Atom *atom = mol.getAtomWithIdx(atomIdx);
atomOrders.push_back(atomIdx);
colors[atomIdx] = GREY_NODE;
// ---------------------
//
// Build the list of possible destinations from here
//
// ---------------------
std::vector<PossibleType> possibles;
possibles.resize(0);
ROMol::OBOND_ITER_PAIR bondsPair = mol.getAtomBonds(atom);
possibles.reserve(bondsPair.second - bondsPair.first);
while (bondsPair.first != bondsPair.second) {
Bond *theBond = mol[*(bondsPair.first)];
bondsPair.first++;
if (bondsInPlay && !(*bondsInPlay)[theBond->getIdx()]) {
continue;
}
if (inBondIdx < 0 ||
theBond->getIdx() != static_cast<unsigned int>(inBondIdx)) {
int otherIdx = theBond->getOtherAtomIdx(atomIdx);
long rank = ranks[otherIdx];
// ---------------------
//
// things are a bit more complicated if we are sitting on a
// ring atom. we would like to traverse first to the
// ring-closure atoms, then to atoms outside the ring first,
// then to atoms in the ring that haven't already been visited
// (non-ring-closure atoms).
//
// Here's how the black magic works:
// - non-ring atom neighbors have their original ranks
// - ring atom neighbors have this added to their ranks:
// (MAX_BONDTYPE - bondOrder)*MAX_NATOMS*MAX_NATOMS
// - ring-closure neighbors lose a factor of:
// (MAX_BONDTYPE+1)*MAX_NATOMS*MAX_NATOMS
//
// This tactic biases us to traverse to non-ring neighbors first,
// original ordering if bond orders are all equal... crafty, neh?
//
// ---------------------
if (!doRandom) {
if (colors[otherIdx] == GREY_NODE) {
rank -= static_cast<int>(MAX_BONDTYPE + 1) * MAX_NATOMS * MAX_NATOMS;
if (!bondSymbols) {
rank += static_cast<int>(MAX_BONDTYPE - theBond->getBondType()) *
MAX_NATOMS;
} else {
const std::string &symb = (*bondSymbols)[theBond->getIdx()];
std::uint32_t hsh = gboost::hash_range(symb.begin(), symb.end());
rank += (hsh % MAX_NATOMS) * MAX_NATOMS;
}
} else if (theBond->getOwningMol().getRingInfo()->numBondRings(
theBond->getIdx())) {
if (!bondSymbols) {
rank += static_cast<int>(MAX_BONDTYPE - theBond->getBondType()) *
MAX_NATOMS * MAX_NATOMS;
} else {
const std::string &symb = (*bondSymbols)[theBond->getIdx()];
std::uint32_t hsh = gboost::hash_range(symb.begin(), symb.end());
rank += (hsh % MAX_NATOMS) * MAX_NATOMS * MAX_NATOMS;
}
}
} else {
// randomize the rank
rank = getRandomGenerator()();
}
// std::cerr << " " << atomIdx << ": " << otherIdx << " " <<
// rank
// << std::endl;
// std::cerr<<"aIdx: "<< atomIdx <<" p: "<<otherIdx<<" Rank:
// "<<ranks[otherIdx] <<" "<<colors[otherIdx]<<"
// "<<theBond->getBondType()<<" "<<rank<<std::endl;
possibles.emplace_back(rank, otherIdx, theBond);
}
}
// ---------------------
//
// Sort on ranks
//
// ---------------------
std::sort(possibles.begin(), possibles.end(), _possibleCompare());
// if (possibles.size())
// std::cerr << " aIdx1: " << atomIdx
// << " first: " << possibles.front().get<0>() << " "
// << possibles.front().get<1>() << std::endl;
// // ---------------------
//
// Now work the children
//
// ---------------------
for (auto &possible : possibles) {
int possibleIdx = possible.get<1>();
Bond *bond = possible.get<2>();
switch (colors[possibleIdx]) {
case WHITE_NODE:
// -----
// we haven't seen this node at all before, traverse
// -----
dfsFindCycles(mol, possibleIdx, bond->getIdx(), colors, ranks,
atomOrders, atomRingClosures, bondsInPlay, bondSymbols,
doRandom);
break;
case GREY_NODE:
// -----
// we've seen this, but haven't finished it (we're finishing a ring)
// -----
atomRingClosures[possibleIdx].push_back(bond->getIdx());
atomRingClosures[atomIdx].push_back(bond->getIdx());
break;
default:
// -----
// this node has been finished. don't do anything.
// -----
break;
}
}
colors[atomIdx] = BLACK_NODE;
} // namespace Canon
void dfsBuildStack(ROMol &mol, int atomIdx, int inBondIdx,
std::vector<AtomColors> &colors, VECT_INT_VECT &cycles,
const UINT_VECT &ranks, UINT_VECT &cyclesAvailable,
MolStack &molStack, UINT_VECT &atomOrders,
UINT_VECT &bondVisitOrders, VECT_INT_VECT &atomRingClosures,
std::vector<INT_LIST> &atomTraversalBondOrder,
const boost::dynamic_bitset<> *bondsInPlay,
const std::vector<std::string> *bondSymbols, bool doRandom) {
#if 0
std::cerr<<"traverse from atom: "<<atomIdx<<" via bond "<<inBondIdx<<" num cycles available: "
<<std::count(cyclesAvailable.begin(),cyclesAvailable.end(),1)<<std::endl;
#endif
Atom *atom = mol.getAtomWithIdx(atomIdx);
INT_LIST directTravList, cycleEndList;
boost::dynamic_bitset<> seenFromHere(mol.getNumAtoms());
seenFromHere.set(atomIdx);
molStack.push_back(MolStackElem(atom));
atomOrders[atom->getIdx()] = rdcast<int>(molStack.size());
colors[atomIdx] = GREY_NODE;
INT_LIST travList;
if (inBondIdx >= 0) {
travList.push_back(inBondIdx);
}
// ---------------------
//
// Add any ring closures
//
// ---------------------
if (atomRingClosures[atomIdx].size()) {
std::vector<unsigned int> ringsClosed;
for (auto bIdx : atomRingClosures[atomIdx]) {
travList.push_back(bIdx);
Bond *bond = mol.getBondWithIdx(bIdx);
seenFromHere.set(bond->getOtherAtomIdx(atomIdx));
unsigned int ringIdx;
if (bond->getPropIfPresent(common_properties::_TraversalRingClosureBond,
ringIdx)) {
// this is end of the ring closure
// we can just pull the ring index from the bond itself:
molStack.push_back(MolStackElem(bond, atomIdx));
bondVisitOrders[bIdx] = molStack.size();
molStack.push_back(MolStackElem(ringIdx));
// don't make the ring digit immediately available again: we don't want
// to have the same
// ring digit opening and closing rings on an atom.
ringsClosed.push_back(ringIdx - 1);
} else {
// this is the beginning of the ring closure, we need to come up with a
// ring index:
auto cAIt =
std::find(cyclesAvailable.begin(), cyclesAvailable.end(), 1);
if (cAIt == cyclesAvailable.end()) {
throw ValueErrorException(
"Too many rings open at once. SMILES cannot be generated.");
}
unsigned int lowestRingIdx = cAIt - cyclesAvailable.begin();
cyclesAvailable[lowestRingIdx] = 0;
++lowestRingIdx;
bond->setProp(common_properties::_TraversalRingClosureBond,
lowestRingIdx);
molStack.push_back(MolStackElem(lowestRingIdx));
}
}
for (auto ringIdx : ringsClosed) {
cyclesAvailable[ringIdx] = 1;
}
}
// ---------------------
//
// Build the list of possible destinations from here
//
// ---------------------
std::vector<PossibleType> possibles;
possibles.resize(0);
ROMol::OBOND_ITER_PAIR bondsPair = mol.getAtomBonds(atom);
possibles.reserve(bondsPair.second - bondsPair.first);
while (bondsPair.first != bondsPair.second) {
Bond *theBond = mol[*(bondsPair.first)];
bondsPair.first++;
if (bondsInPlay && !(*bondsInPlay)[theBond->getIdx()]) {
continue;
}
if (inBondIdx < 0 ||
theBond->getIdx() != static_cast<unsigned int>(inBondIdx)) {
int otherIdx = theBond->getOtherAtomIdx(atomIdx);
// ---------------------
//
// This time we skip the ring-closure atoms (we did them
// above); we want to traverse first to atoms outside the ring
// then to atoms in the ring that haven't already been visited
// (non-ring-closure atoms).
//
// otherwise it's the same ranking logic as above
//
// ---------------------
if (colors[otherIdx] != WHITE_NODE || seenFromHere[otherIdx]) {
// ring closure or finished atom... skip it.
continue;
}
unsigned long rank = ranks[otherIdx];
if (!doRandom) {
if (theBond->getOwningMol().getRingInfo()->numBondRings(
theBond->getIdx())) {
if (!bondSymbols) {
rank += static_cast<int>(MAX_BONDTYPE - theBond->getBondType()) *
MAX_NATOMS * MAX_NATOMS;
} else {
const std::string &symb = (*bondSymbols)[theBond->getIdx()];
std::uint32_t hsh = gboost::hash_range(symb.begin(), symb.end());
rank += (hsh % MAX_NATOMS) * MAX_NATOMS * MAX_NATOMS;
}
}
} else {
// randomize the rank
rank = getRandomGenerator()();
}
possibles.emplace_back(rank, otherIdx, theBond);
}
}
// ---------------------
//
// Sort on ranks
//
// ---------------------
std::sort(possibles.begin(), possibles.end(), _possibleCompare());
// if (possibles.size())
// std::cerr << " aIdx2: " << atomIdx
// << " first: " << possibles.front().get<0>() << " "
// << possibles.front().get<1>() << std::endl;
// ---------------------
//
// Now work the children
//
// ---------------------
for (auto possiblesIt = possibles.begin(); possiblesIt != possibles.end();
possiblesIt++) {
int possibleIdx = possiblesIt->get<1>();
if (colors[possibleIdx] != WHITE_NODE) {
// we're either done or it's a ring-closure, which we already processed...
// this test isn't strictly required, because we only added WHITE notes to
// the possibles list, but it seems logical to document it
continue;
}
Bond *bond = possiblesIt->get<2>();
Atom *otherAtom = mol.getAtomWithIdx(possibleIdx);
// ww might have some residual data from earlier calls, clean that up:
otherAtom->clearProp(common_properties::_TraversalBondIndexOrder);
travList.push_back(bond->getIdx());
if (possiblesIt + 1 != possibles.end()) {
// we're branching
molStack.push_back(
MolStackElem("(", rdcast<int>(possiblesIt - possibles.begin())));
}
molStack.push_back(MolStackElem(bond, atomIdx));
bondVisitOrders[bond->getIdx()] = molStack.size();
dfsBuildStack(mol, possibleIdx, bond->getIdx(), colors, cycles, ranks,
cyclesAvailable, molStack, atomOrders, bondVisitOrders,
atomRingClosures, atomTraversalBondOrder, bondsInPlay,
bondSymbols, doRandom);
if (possiblesIt + 1 != possibles.end()) {
molStack.push_back(
MolStackElem(")", rdcast<int>(possiblesIt - possibles.begin())));
}
}
atomTraversalBondOrder[atom->getIdx()] = travList;
colors[atomIdx] = BLACK_NODE;
}
void canonicalDFSTraversal(ROMol &mol, int atomIdx, int inBondIdx,
std::vector<AtomColors> &colors,
VECT_INT_VECT &cycles, const UINT_VECT &ranks,
UINT_VECT &cyclesAvailable, MolStack &molStack,
UINT_VECT &atomOrders, UINT_VECT &bondVisitOrders,
VECT_INT_VECT &atomRingClosures,
std::vector<INT_LIST> &atomTraversalBondOrder,
const boost::dynamic_bitset<> *bondsInPlay,
const std::vector<std::string> *bondSymbols,
bool doRandom) {
PRECONDITION(colors.size() >= mol.getNumAtoms(), "vector too small");
PRECONDITION(ranks.size() >= mol.getNumAtoms(), "vector too small");
PRECONDITION(atomOrders.size() >= mol.getNumAtoms(), "vector too small");
PRECONDITION(bondVisitOrders.size() >= mol.getNumBonds(), "vector too small");
PRECONDITION(atomRingClosures.size() >= mol.getNumAtoms(),
"vector too small");
PRECONDITION(atomTraversalBondOrder.size() >= mol.getNumAtoms(),
"vector too small");
PRECONDITION(!bondsInPlay || bondsInPlay->size() >= mol.getNumBonds(),
"bondsInPlay too small");
PRECONDITION(!bondSymbols || bondSymbols->size() >= mol.getNumBonds(),
"bondSymbols too small");
std::vector<AtomColors> tcolors;
tcolors.resize(colors.size());
std::copy(colors.begin(), colors.end(), tcolors.begin());
dfsFindCycles(mol, atomIdx, inBondIdx, tcolors, ranks, atomOrders,
atomRingClosures, bondsInPlay, bondSymbols, doRandom);
dfsBuildStack(mol, atomIdx, inBondIdx, colors, cycles, ranks, cyclesAvailable,
molStack, atomOrders, bondVisitOrders, atomRingClosures,
atomTraversalBondOrder, bondsInPlay, bondSymbols, doRandom);
}
bool canHaveDirection(const Bond *bond) {
PRECONDITION(bond, "bad bond");
Bond::BondType bondType = bond->getBondType();
return (bondType == Bond::SINGLE || bondType == Bond::AROMATIC);
}
void clearBondDirs(ROMol &mol, Bond *refBond, const Atom *fromAtom,
UINT_VECT &bondDirCounts, UINT_VECT &atomDirCounts,
const UINT_VECT &bondVisitOrders) {
RDUNUSED_PARAM(bondVisitOrders);
PRECONDITION(bondDirCounts.size() >= mol.getNumBonds(), "bad dirCount size");
PRECONDITION(refBond, "bad bond");
PRECONDITION(&refBond->getOwningMol() == &mol, "bad bond");
PRECONDITION(fromAtom, "bad atom");
PRECONDITION(&fromAtom->getOwningMol() == &mol, "bad bond");
#if 0
std::copy(bondDirCounts.begin(),bondDirCounts.end(),std::ostream_iterator<int>(std::cerr,", "));
std::cerr<<"\n";
std::copy(atomDirCounts.begin(),atomDirCounts.end(),std::ostream_iterator<int>(std::cerr,", "));
std::cerr<<"\n";
std::cerr<<"cBD: bond: "<<refBond->getIdx()<<" atom: "<<fromAtom->getIdx()<<": ";
#endif
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(fromAtom);
bool nbrPossible = false, adjusted = false;
while (beg != end) {
Bond *oBond = mol[*beg];
// std::cerr<<" >>"<<oBond->getIdx()<<" "<<canHaveDirection(oBond)<<"
// "<<bondDirCounts[oBond->getIdx()]<<"-"<<bondDirCounts[refBond->getIdx()]<<"
// "<<atomDirCounts[oBond->getBeginAtomIdx()]<<"-"<<atomDirCounts[oBond->getEndAtomIdx()]<<std::endl;
if (oBond != refBond && canHaveDirection(oBond)) {
nbrPossible = true;
if ((bondDirCounts[oBond->getIdx()] >=
bondDirCounts[refBond->getIdx()]) &&
atomDirCounts[oBond->getBeginAtomIdx()] != 1 &&
atomDirCounts[oBond->getEndAtomIdx()] != 1) {
adjusted = true;
bondDirCounts[oBond->getIdx()] -= 1;
if (!bondDirCounts[oBond->getIdx()]) {
// no one is setting the direction here:
oBond->setBondDir(Bond::NONE);
atomDirCounts[oBond->getBeginAtomIdx()] -= 1;
atomDirCounts[oBond->getEndAtomIdx()] -= 1;
// std::cerr<<"ob:"<<oBond->getIdx()<<" ";
}
}
}
beg++;
}
if (nbrPossible && !adjusted &&
atomDirCounts[refBond->getBeginAtomIdx()] != 1 &&
atomDirCounts[refBond->getEndAtomIdx()] != 1) {
// we found a neighbor that could have directionality set,
// but it had a lower bondDirCount than us, so we must
// need to be adjusted:
bondDirCounts[refBond->getIdx()] -= 1;
if (!bondDirCounts[refBond->getIdx()]) {
refBond->setBondDir(Bond::NONE);
atomDirCounts[refBond->getBeginAtomIdx()] -= 1;
atomDirCounts[refBond->getEndAtomIdx()] -= 1;
// std::cerr<<"rb:"<<refBond->getIdx()<<" ";
}
}
// std::cerr<<std::endl;
}
void removeRedundantBondDirSpecs(ROMol &mol, MolStack &molStack,
UINT_VECT &bondDirCounts,
UINT_VECT &atomDirCounts,
const UINT_VECT &bondVisitOrders) {
PRECONDITION(bondDirCounts.size() >= mol.getNumBonds(), "bad dirCount size");
#if 0
std::cerr<<"rRBDS: ";
mol.debugMol(std::cerr);
std::copy(bondDirCounts.begin(),bondDirCounts.end(),std::ostream_iterator<int>(std::cerr,", "));
std::cerr<<"\n";
#endif
// find bonds that have directions indicated that are redundant:
for (auto &msI : molStack) {
if (msI.type == MOL_STACK_BOND) {
Bond *tBond = msI.obj.bond;
const Atom *canonBeginAtom = mol.getAtomWithIdx(msI.number);
const Atom *canonEndAtom =
mol.getAtomWithIdx(tBond->getOtherAtomIdx(msI.number));
if (canHaveDirection(tBond) && bondDirCounts[tBond->getIdx()] >= 1) {
// start by finding the double bond that sets tBond's direction:
const Atom *dblBondAtom = nullptr;
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(canonBeginAtom);
while (beg != end) {
if (mol[*beg] != tBond && mol[*beg]->getBondType() == Bond::DOUBLE &&
mol[*beg]->getStereo() > Bond::STEREOANY) {
dblBondAtom =
canonBeginAtom; // tBond->getOtherAtom(canonBeginAtom);
break;
}
beg++;
}
if (dblBondAtom != nullptr) {
clearBondDirs(mol, tBond, dblBondAtom, bondDirCounts, atomDirCounts,
bondVisitOrders);
}
dblBondAtom = nullptr;
boost::tie(beg, end) = mol.getAtomBonds(canonEndAtom);
while (beg != end) {
if (mol[*beg] != tBond && mol[*beg]->getBondType() == Bond::DOUBLE &&
mol[*beg]->getStereo() > Bond::STEREOANY) {
dblBondAtom = canonEndAtom; // tBond->getOtherAtom(canonEndAtom);
break;
}
beg++;
}
if (dblBondAtom != nullptr) {
clearBondDirs(mol, tBond, dblBondAtom, bondDirCounts, atomDirCounts,
bondVisitOrders);
}
} else if (tBond->getBondDir() != Bond::NONE) {
// we aren't supposed to have a direction set, but we do:
tBond->setBondDir(Bond::NONE);
}
}
}
}
void canonicalizeFragment(ROMol &mol, int atomIdx,
std::vector<AtomColors> &colors,
const UINT_VECT &ranks, MolStack &molStack,
const boost::dynamic_bitset<> *bondsInPlay,
const std::vector<std::string> *bondSymbols,
bool doIsomericSmiles, bool doRandom) {
PRECONDITION(colors.size() >= mol.getNumAtoms(), "vector too small");
PRECONDITION(ranks.size() >= mol.getNumAtoms(), "vector too small");
PRECONDITION(!bondsInPlay || bondsInPlay->size() >= mol.getNumBonds(),
"bondsInPlay too small");
PRECONDITION(!bondSymbols || bondSymbols->size() >= mol.getNumBonds(),
"bondSymbols too small");
unsigned int nAtoms = mol.getNumAtoms();
UINT_VECT atomVisitOrders(nAtoms, 0);
UINT_VECT bondVisitOrders(mol.getNumBonds(), 0);
UINT_VECT bondDirCounts(mol.getNumBonds(), 0);
UINT_VECT atomDirCounts(nAtoms, 0);
UINT_VECT cyclesAvailable(MAX_CYCLES, 1);
VECT_INT_VECT cycles(nAtoms);
boost::dynamic_bitset<> ringStereoChemAdjusted(nAtoms);
// make sure that we've done the stereo perception:
if (!mol.hasProp(common_properties::_StereochemDone)) {
MolOps::assignStereochemistry(mol, false);
}
// we need ring information; make sure findSSSR has been called before
// if not call now
if (!mol.getRingInfo()->isInitialized()) {
MolOps::findSSSR(mol);
}
mol.getAtomWithIdx(atomIdx)->setProp(common_properties::_TraversalStartPoint,
true);
VECT_INT_VECT atomRingClosures(nAtoms);
std::vector<INT_LIST> atomTraversalBondOrder(nAtoms);
Canon::canonicalDFSTraversal(
mol, atomIdx, -1, colors, cycles, ranks, cyclesAvailable, molStack,
atomVisitOrders, bondVisitOrders, atomRingClosures,
atomTraversalBondOrder, bondsInPlay, bondSymbols, doRandom);
CHECK_INVARIANT(!molStack.empty(), "Empty stack.");
CHECK_INVARIANT(molStack.begin()->type == MOL_STACK_ATOM,
"Corrupted stack. First element should be an atom.");
// collect some information about traversal order on chiral atoms
boost::dynamic_bitset<> numSwapsChiralAtoms(nAtoms);
if (doIsomericSmiles) {
for (const auto atom : mol.atoms()) {
if (atom->getChiralTag() != Atom::CHI_UNSPECIFIED) {
// check if all of this atom's bonds are in play
for (const auto &bndItr :
boost::make_iterator_range(mol.getAtomBonds(atom))) {
if (bondsInPlay && !(*bondsInPlay)[mol[bndItr]->getIdx()]) {
atom->setProp(common_properties::_brokenChirality, true);
break;
}
}
if (atom->hasProp(common_properties::_brokenChirality)) {
continue;
}
const INT_LIST &trueOrder = atomTraversalBondOrder[atom->getIdx()];
// Check if the atom can be chiral, and if chirality needs inversion
if (trueOrder.size() >= 3) {
int nSwaps;
// We have to make sure that trueOrder contains all the
// bonds, even if they won't be written to the SMARTS
if (trueOrder.size() < atom->getDegree()) {
INT_LIST tOrder = trueOrder;
for (const auto &bndItr :
boost::make_iterator_range(mol.getAtomBonds(atom))) {
int bndIdx = mol[bndItr]->getIdx();
if (std::find(trueOrder.begin(), trueOrder.end(), bndIdx) ==
trueOrder.end()) {
tOrder.push_back(bndIdx);
break;
}
}
nSwaps = atom->getPerturbationOrder(tOrder);
} else {
nSwaps = atom->getPerturbationOrder(trueOrder);
}
if (chiralAtomNeedsTagInversion(
mol, atom,
molStack.begin()->obj.atom->getIdx() == atom->getIdx(),
atomRingClosures[atom->getIdx()].size())) {
// This is a special case. Here's an example:
// Our internal representation of a chiral center is equivalent
// to:
// [C@](F)(O)(C)[H]
// we'll be dumping it without the H, which entails a
// reordering:
// [C@@H](F)(O)C
++nSwaps;
}
if (nSwaps % 2) {
numSwapsChiralAtoms.set(atom->getIdx());
}
}
}
}
}
// remove the current directions on single bonds around double bonds:
for (auto bond : mol.bonds()) {
Bond::BondDir dir = bond->getBondDir();
if (dir == Bond::ENDDOWNRIGHT || dir == Bond::ENDUPRIGHT) {
bond->setBondDir(Bond::NONE);
}
}
#if 0
std::cerr<<"<11111111"<<std::endl;
std::cerr<<"----------------------------------------->"<<std::endl;
mol.debugMol(std::cerr);
#endif
// std::cerr<<"----->\ntraversal stack:"<<std::endl;
// traverse the stack and canonicalize double bonds and atoms with (ring)
// stereochemistry
for (auto &msI : molStack) {
#if 0
if(msI->type == MOL_STACK_ATOM) std::cerr<<" atom: "<<msI->obj.atom->getIdx()<<std::endl;
else if(msI->type == MOL_STACK_BOND) std::cerr<<" bond: "<<msI->obj.bond->getIdx()<<" "<<msI->number<<" "<<msI->obj.bond->getBeginAtomIdx()<<"-"<<msI->obj.bond->getEndAtomIdx()<<" order: "<<msI->obj.bond->getBondType()<<std::endl;
else if(msI->type == MOL_STACK_RING) std::cerr<<" ring: "<<msI->number<<std::endl;
else if(msI->type == MOL_STACK_BRANCH_OPEN) std::cerr<<" branch open"<<std::endl;
else if(msI->type == MOL_STACK_BRANCH_CLOSE) std::cerr<<" branch close"<<std::endl;
#endif
if (msI.type == MOL_STACK_BOND &&
msI.obj.bond->getBondType() == Bond::DOUBLE &&
msI.obj.bond->getStereo() > Bond::STEREOANY) {
if (msI.obj.bond->getStereoAtoms().size() >= 2) {
Canon::canonicalizeDoubleBond(msI.obj.bond, bondVisitOrders,
atomVisitOrders, bondDirCounts,
atomDirCounts, molStack);
} else {
// bad stereo spec:
msI.obj.bond->setStereo(Bond::STEREONONE);
}
}
if (doIsomericSmiles) {
if (msI.type == MOL_STACK_ATOM &&
msI.obj.atom->getChiralTag() != Atom::CHI_UNSPECIFIED &&
!msI.obj.atom->hasProp(common_properties::_brokenChirality)) {
if (msI.obj.atom->hasProp(common_properties::_ringStereoAtoms)) {
if (!ringStereoChemAdjusted[msI.obj.atom->getIdx()]) {
msI.obj.atom->setChiralTag(Atom::CHI_TETRAHEDRAL_CCW);
ringStereoChemAdjusted.set(msI.obj.atom->getIdx());
}
const INT_VECT &ringStereoAtoms = msI.obj.atom->getProp<INT_VECT>(
common_properties::_ringStereoAtoms);
for (auto nbrV : ringStereoAtoms) {
int nbrIdx = abs(nbrV) - 1;
// Adjust the chirality flag of the ring stereo atoms according to
// the first one
if (!ringStereoChemAdjusted[nbrIdx] &&
atomVisitOrders[nbrIdx] >
atomVisitOrders[msI.obj.atom->getIdx()]) {
mol.getAtomWithIdx(nbrIdx)->setChiralTag(
msI.obj.atom->getChiralTag());
if (nbrV < 0) {
mol.getAtomWithIdx(nbrIdx)->invertChirality();
}
// Odd number of swaps for first chiral ring atom --> needs to be
// swapped but we want to retain chirality
if (numSwapsChiralAtoms[msI.obj.atom->getIdx()]) {
// Odd number of swaps for chiral ring neighbor --> needs to be
// swapped but we want to retain chirality
if (!numSwapsChiralAtoms[nbrIdx]) {
mol.getAtomWithIdx(nbrIdx)->invertChirality();
}
}
// Even number of swaps for first chiral ring atom --> don't need
// to be swapped
else {
// Odd number of swaps for chiral ring neighbor --> needs to be
// swapped
if (numSwapsChiralAtoms[nbrIdx]) {
mol.getAtomWithIdx(nbrIdx)->invertChirality();
}
}
ringStereoChemAdjusted.set(nbrIdx);
}
}
} else {
if (msI.obj.atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW) {
if ((numSwapsChiralAtoms[msI.obj.atom->getIdx()])) {
msI.obj.atom->invertChirality();
}
} else if (msI.obj.atom->getChiralTag() ==
Atom::CHI_TETRAHEDRAL_CCW) {
if ((numSwapsChiralAtoms[msI.obj.atom->getIdx()])) {
msI.obj.atom->invertChirality();
}
}
}
}
}
}
#if 0
std::cerr<<"<-----"<<std::endl;
std::cerr<<"----------------------------------------->"<<std::endl;
mol.debugMol(std::cerr);
#endif
Canon::removeRedundantBondDirSpecs(mol, molStack, bondDirCounts,
atomDirCounts, bondVisitOrders);
#if 0
std::cerr<<"----------------------------------------->"<<std::endl;
mol.debugMol(std::cerr);
std::cerr<<"----------------------------------------->"<<std::endl;
#endif
}
} // namespace Canon
} // namespace RDKit
Reference in New Issue View Git Blame Copy Permalink