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
2d257b11944303ee161aea564e9c3bce276df6cc
rdkit/Code/GraphMol/Canon.cpp
Hussein Faara f35e7e6414 remove no-op macros and dead code (pt 1) (#8012)
2024-11-19 04:31:56 +01:00

1367 lines
52 KiB
C++
Raw Blame History

//
// Copyright (C) 2001-2021 Greg Landrum and other RDKit contributors
//
// @@ 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/Chirality.h>
#include <GraphMol/new_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 details {
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) {
auto cDescr = (*cIt)->getDescription();
if (cDescr == "AtomHCount") {
return !(*cIt)->getNegation() &&
((ATOM_EQUALS_QUERY *)(*cIt).get())->getVal() == 1;
} else if (cDescr == "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;
}
} // namespace details
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) ||
(!details::atomHasFourthValence(atom) && numClosures == 1 &&
!details::isUnsaturated(atom, mol)));
}
auto _possibleCompare = [](const PossibleType &arg1, const PossibleType &arg2) {
return (std::get<0>(arg1) < std::get<0>(arg2));
};
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() ||
isDative(*bond),
"bad bond type");
switch (bond->getBondDir()) {
case Bond::ENDUPRIGHT:
bond->setBondDir(Bond::ENDDOWNRIGHT);
break;
case Bond::ENDDOWNRIGHT:
bond->setBondDir(Bond::ENDUPRIGHT);
break;
default:
break;
}
}
namespace {
bool 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);
}
} // namespace
// 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;
auto findNeighborBonds = [&mol, &dblBond, &bondDirCounts, &bondVisitOrders,
&firstVisitOrder](
auto atom, auto &firstNeighborBond,
auto &secondNeighborBond, auto &dirSet) {
for (const auto bond : mol.atomBonds(atom)) {
if (bond == dblBond || !canSetDoubleBondStereo(*bond)) {
continue;
}
auto bondIdx = bond->getIdx();
if (bondDirCounts[bondIdx] > 0) {
dirSet = true;
}
if (!firstNeighborBond || bondVisitOrders[bondIdx] < firstVisitOrder) {
if (firstNeighborBond) {
secondNeighborBond = firstNeighborBond;
}
firstNeighborBond = bond;
firstVisitOrder = bondVisitOrders[bondIdx];
} else {
secondNeighborBond = bond;
}
}
};
findNeighborBonds(atom1, firstFromAtom1, secondFromAtom1, dir1Set);
firstVisitOrder = mol.getNumBonds() + 1;
findNeighborBonds(atom2, firstFromAtom2, secondFromAtom2, dir2Set);
// Make sure we found everything we need to find.
// This really shouldn't be a problem, but molecules can end up in odd
// states; for example, allenes can end up here. Instead of checking for them
// explicitly, exit early in any such possible state.
if (!firstFromAtom1 || !firstFromAtom2) {
return;
}
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;
} 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
auto flipBondDir = [](Bond::BondDir bondDir) {
return (bondDir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT
: Bond::ENDUPRIGHT;
};
// 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 = flipBondDir(atom1Dir);
}
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 isFlipped = false;
if (atom1->getDegree() == 3 && // atom1ControllingBond == firstFromAtom1 &&
std::find(stereoAtoms.begin(), stereoAtoms.end(),
static_cast<int>(atom1ControllingBond->getOtherAtomIdx(
atom1->getIdx()))) == stereoAtoms.end()) {
isFlipped = true;
atom2Dir = flipBondDir(atom2Dir);
}
// 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 = flipBondDir(atom2Dir);
}
if (!isFlipped && isClosingRingBond(dblBond)) {
atom2Dir = flipBondDir(atom2Dir);
}
// 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 = flipBondDir(atom2Dir);
}
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 = flipBondDir(atom1Dir);
}
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 = flipBondDir(atom1Dir);
}
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 && secondFromAtom1) {
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 = flipBondDir(firstFromAtom1->getBondDir());
secondFromAtom1->setBondDir(otherDir);
} else {
secondFromAtom1->setBondDir(firstFromAtom1->getBondDir());
}
}
bondDirCounts[secondFromAtom1->getIdx()] += 1;
atomDirCounts[atom1->getIdx()] += 1;
}
if (atom2->getDegree() == 3 && secondFromAtom2) {
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 = flipBondDir(firstFromAtom2->getBondDir());
} 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;
auto dir = firstFromAtom2->getBondDir();
if (isClosingRingBond(otherAtom3Bond)) {
dir = flipBondDir(dir);
}
otherAtom3Bond->setBondDir(dir);
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()std:std::get<0>() << " "
// << possibles.front()std:std::get<1>() << std::endl;
// // ---------------------
//
// Now work the children
//
// ---------------------
for (auto &possible : possibles) {
int possibleIdx = std::get<1>(possible);
Bond *bond = std::get<2>(possible);
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) {
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()std:std::get<0>() << " "
// << possibles.front()std:std::get<1>() << std::endl;
// ---------------------
//
// Now work the children
//
// ---------------------
for (auto possiblesIt = possibles.begin(); possiblesIt != possibles.end();
possiblesIt++) {
int possibleIdx = std::get<1>(*possiblesIt);
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 = std::get<2>(*possiblesIt);
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);
}
void clearBondDirs(ROMol &mol, Bond *refBond, const Atom *fromAtom,
UINT_VECT &bondDirCounts, UINT_VECT &atomDirCounts,
const UINT_VECT &) {
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");
ROMol::OEDGE_ITER beg, end;
boost::tie(beg, end) = mol.getAtomBonds(fromAtom);
bool nbrPossible = false, adjusted = false;
while (beg != end) {
Bond *oBond = mol[*beg];
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");
// 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);
}
}
}
}
// insert (-1) for hydrogens or missing ligands, where these are placed
// depends on if it is the first atom or not
static void insertImplicitNbors(INT_LIST &bonds, const Atom::ChiralType tag,
const bool firstAtom) {
unsigned int ref_max = Chirality::getMaxNbors(tag);
if (bonds.size() < ref_max) {
if (firstAtom) {
bonds.insert(bonds.begin(), ref_max - bonds.size(), -1);
} else {
bonds.insert(++bonds.begin(), ref_max - bonds.size(), -1);
}
}
}
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,
bool doChiralInversions) {
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
// NOTE: if called from the SMARTS code, the ring info will be set to SSSR,
// but no ring infor in actually set
if (!mol.getRingInfo()->isSymmSssr()) {
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);
std::vector<int> atomPermutationIndices(nAtoms, 0);
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 bnd : mol.atomBonds(atom)) {
if (bondsInPlay && !(*bondsInPlay)[bnd->getIdx()]) {
atom->setProp(common_properties::_brokenChirality, true);
break;
}
}
if (atom->hasProp(common_properties::_brokenChirality)) {
continue;
}
// Check if the atom can be chiral, and if chirality needs inversion
const INT_LIST &trueOrder = atomTraversalBondOrder[atom->getIdx()];
// Extra check needed if/when @AL1/@AL2 supported
if (trueOrder.size() >= 3 || Chirality::hasNonTetrahedralStereo(atom)) {
int nSwaps = 0;
int perm = 0;
if (Chirality::hasNonTetrahedralStereo(atom)) {
atom->getPropIfPresent(common_properties::_chiralPermutation, perm);
}
const unsigned int firstIdx = molStack.begin()->obj.atom->getIdx();
const bool firstInPart = atom->getIdx() == firstIdx;
// We have to make sure that trueOrder contains all the
// bonds, even if they won't be written to the SMILES
if (trueOrder.size() < atom->getDegree()) {
INT_LIST tOrder = trueOrder;
for (const auto bnd : mol.atomBonds(atom)) {
int bndIdx = bnd->getIdx();
if (std::find(trueOrder.begin(), trueOrder.end(), bndIdx) ==
trueOrder.end()) {
tOrder.push_back(bndIdx);
break;
}
}
if (!perm) {
nSwaps = atom->getPerturbationOrder(tOrder);
} else {
insertImplicitNbors(tOrder, atom->getChiralTag(), firstInPart);
perm = Chirality::getChiralPermutation(atom, tOrder);
}
} else {
if (!perm) {
nSwaps = atom->getPerturbationOrder(trueOrder);
} else {
INT_LIST tOrder = trueOrder;
insertImplicitNbors(tOrder, atom->getChiralTag(), firstInPart);
perm = Chirality::getChiralPermutation(atom, tOrder);
}
}
// in future this should be moved up and simplified, there should not
// be an option to not do chiral inversions
if (doChiralInversions &&
chiralAtomNeedsTagInversion(
mol, atom, firstInPart,
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());
}
atomPermutationIndices[atom->getIdx()] = perm;
}
}
}
}
// 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);
}
}
// traverse the stack and canonicalize double bonds and atoms with (ring)
// stereochemistry
for (auto &msI : molStack) {
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)) {
// FIX: handle stereogroups here too
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 (size_t sgidx;
msI.obj.atom->getPropIfPresent("_stereoGroup", sgidx) &&
mol.getStereoGroups().size() > sgidx) {
// make sure that the reference atom in the stereogroup is CCW
auto &sg = mol.getStereoGroups()[sgidx];
bool swapIt =
msI.obj.atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW;
if (swapIt) {
msI.obj.atom->invertChirality();
}
if (swapIt || numSwapsChiralAtoms[msI.obj.atom->getIdx()]) {
for (auto at : sg.getAtoms()) {
if (at == msI.obj.atom) {
continue;
}
at->invertChirality();
}
}
} else {
if (msI.obj.atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW ||
msI.obj.atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CCW) {
if ((numSwapsChiralAtoms[msI.obj.atom->getIdx()])) {
msI.obj.atom->invertChirality();
}
} else if (atomPermutationIndices[msI.obj.atom->getIdx()]) {
msI.obj.atom->setProp(
common_properties::_chiralPermutation,
atomPermutationIndices[msI.obj.atom->getIdx()]);
}
}
}
}
}
Canon::removeRedundantBondDirSpecs(mol, molStack, bondDirCounts,
atomDirCounts, bondVisitOrders);
}
void canonicalizeEnhancedStereo(ROMol &mol,
const std::vector<unsigned int> *atomRanks) {
const auto &sgs = mol.getStereoGroups();
if (sgs.empty()) {
return;
}
std::vector<unsigned int> lranks;
if (!atomRanks) {
bool breakTies = true;
rankMolAtoms(mol, lranks, breakTies);
atomRanks = &lranks;
}
// one thing that makes this all easier is that the stereogroups are
// independent of each other
std::vector<StereoGroup> newSgs;
for (auto &sg : sgs) {
// we don't do anything to ABS groups
if (sg.getGroupType() == StereoGroupType::STEREO_ABSOLUTE) {
newSgs.push_back(sg);
continue;
}
// sort the atoms by rank:
auto getAtomRank = [&atomRanks](const Atom *at1, const Atom *at2) {
return atomRanks->at(at1->getIdx()) < atomRanks->at(at2->getIdx());
};
auto sgAtoms = sg.getAtoms();
std::sort(sgAtoms.begin(), sgAtoms.end(), getAtomRank);
// sort the bonds by atom rank:
auto getBondRank = [&atomRanks](const Bond *bd1, const Bond *bd2) {
unsigned int bd1at1 = atomRanks->at(bd1->getBeginAtomIdx());
unsigned int bd1at2 = atomRanks->at(bd1->getEndAtomIdx());
unsigned int bd2at1 = atomRanks->at(bd2->getBeginAtomIdx());
unsigned int bd2at2 = atomRanks->at(bd2->getEndAtomIdx());
if (bd1at1 < bd1at2) {
std::swap(bd1at1, bd1at2);
}
if (bd2at1 < bd2at2) {
std::swap(bd2at1, bd2at2);
}
if (bd1at1 != bd2at1) {
return bd1at1 < bd2at1;
}
return bd1at2 < bd2at2;
};
auto sgBonds = sg.getBonds();
std::sort(sgBonds.begin(), sgBonds.end(), getBondRank);
// find the reference (lowest-ranked) atom (or lowest-ranked bond)
Atom::ChiralType foundRefState = Atom::ChiralType::CHI_TETRAHEDRAL_CCW;
if (sgAtoms.size() > 0) {
foundRefState = sgAtoms.front()->getChiralTag();
} else if (sgBonds.size() > 0) {
if (sgBonds.front()->getStereo() == Bond::BondStereo::STEREOATROPCCW) {
foundRefState =
Atom::ChiralType::CHI_TETRAHEDRAL_CCW; // convert atropisomer CCW
// to atom CCW
} else {
foundRefState =
Atom::ChiralType::CHI_TETRAHEDRAL_CW; // convert atropisomer CW
// to atom CW
}
}
// we will use CCW as the "canonical" state for chirality, so if the
// referenceAtom is already CCW then we don't need to do anything more
// with this stereogroup
auto refState = Atom::ChiralType::CHI_TETRAHEDRAL_CCW;
if (foundRefState != refState) {
// we need to flip everyone... so loop over the other atoms and bonds
// and flip them all:
for (auto atom : sgAtoms) {
atom->invertChirality();
}
for (auto bond : sgBonds) {
bond->invertChirality();
}
}
newSgs.emplace_back(
StereoGroup(sg.getGroupType(), std::move(sgAtoms), std::move(sgBonds)));
// note that we do not forward the Group Ids: this is intentional, so that
// the Ids are reassigned based on the canonicalized order.
if (sgAtoms.size() > 0) {
sgAtoms.front()->setProp("_stereoGroup", newSgs.size() - 1, true);
}
}
mol.setStereoGroups(newSgs);
}
void addSingleAbsGroup(ROMol &mol) {
// all chiral centers are added to an abs group
// if there are not chiral centers, no group is added
std::vector<StereoGroup> sgs;
std::vector<Atom *> chiralAtoms;
std::vector<Bond *> chiralBonds;
for (auto &atom : mol.atoms()) {
if (atom->getChiralTag() == Atom::ChiralType::CHI_TETRAHEDRAL_CCW ||
atom->getChiralTag() == Atom::ChiralType::CHI_TETRAHEDRAL_CW) {
chiralAtoms.push_back(atom);
}
}
for (auto &bond : mol.bonds()) {
if (bond->getStereo() == Bond::BondStereo::STEREOATROPCW ||
bond->getStereo() == Bond::BondStereo::STEREOATROPCCW) {
chiralBonds.push_back(bond);
}
}
if (!chiralAtoms.empty() || !chiralBonds.empty()) {
sgs.emplace_back(StereoGroupType::STEREO_ABSOLUTE, chiralAtoms,
chiralBonds);
}
mol.setStereoGroups(sgs); // could be empty, or have one abs group
}
void clearStereoGroups(ROMol &mol) {
// all chiral centers are added to an abs group
// if there are not chiral centers, no group is added
std::vector<StereoGroup> sgs;
mol.setStereoGroups(sgs);
}
}; // namespace Canon
} // namespace RDKit
Reference in New Issue View Git Blame Copy Permalink