mirror of
https://github.com/rdkit/rdkit.git
synced 2026-06-03 21:44:30 +08:00
ecf5a67acc
* Fixes #2437 getting the canonical atom ranking no longer results in molecules have a RingInfo structure that's been initialized but contains nothing. * update expected results for the MMPA tests
642 lines
22 KiB
C++
642 lines
22 KiB
C++
//
|
|
// Copyright (C) 2014 Greg Landrum
|
|
// Adapted from pseudo-code from Roger Sayle
|
|
//
|
|
// @@ 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 "new_canon.h"
|
|
#include <GraphMol/RDKitBase.h>
|
|
#include <cstdint>
|
|
#include <boost/foreach.hpp>
|
|
#include <cstring>
|
|
#include <iostream>
|
|
#include <cassert>
|
|
|
|
namespace RDKit {
|
|
namespace Canon {
|
|
void CreateSinglePartition(unsigned int nAtoms, int *order, int *count,
|
|
canon_atom *atoms) {
|
|
for (unsigned int i = 0; i < nAtoms; i++) {
|
|
atoms[i].index = 0;
|
|
order[i] = i;
|
|
count[i] = 0;
|
|
}
|
|
count[0] = nAtoms;
|
|
}
|
|
|
|
void ActivatePartitions(unsigned int nAtoms, int *order, int *count,
|
|
int &activeset, int *next, int *changed) {
|
|
unsigned int i, j;
|
|
activeset = -1;
|
|
for (i = 0; i < nAtoms; i++) next[i] = -2;
|
|
|
|
i = 0;
|
|
do {
|
|
j = order[i];
|
|
if (count[j] > 1) {
|
|
next[j] = activeset;
|
|
activeset = j;
|
|
i += count[j];
|
|
} else
|
|
i++;
|
|
} while (i < nAtoms);
|
|
|
|
for (i = 0; i < nAtoms; i++) {
|
|
j = order[i];
|
|
int flag = 1;
|
|
//#define SKIP_NODE_CHANGED_OPTIMIZATION 0
|
|
//#ifndef SKIP_NODE_CHANGED_OPTIMIZATION
|
|
// if(count[j]){
|
|
// std::cout << "j " << j << std::endl;
|
|
// flag=(next[j]!=-2);
|
|
// }
|
|
//#endif
|
|
changed[j] = flag;
|
|
}
|
|
}
|
|
|
|
void compareRingAtomsConcerningNumNeighbors(Canon::canon_atom *atoms,
|
|
unsigned int nAtoms,
|
|
const ROMol &mol) {
|
|
RingInfo *ringInfo = mol.getRingInfo();
|
|
for (unsigned idx = 0; idx < nAtoms; ++idx) {
|
|
const Canon::canon_atom &a = atoms[idx];
|
|
if (!ringInfo->isInitialized() || ringInfo->numAtomRings(a.atom->getIdx()) < 1) {
|
|
continue;
|
|
}
|
|
std::deque<int> neighbors;
|
|
neighbors.push_back(idx);
|
|
unsigned currentRNIdx = 0;
|
|
atoms[idx].neighborNum.reserve(1000);
|
|
atoms[idx].revistedNeighbors.assign(1000, 0);
|
|
char *visited = (char *)malloc(nAtoms * sizeof(char));
|
|
memset(visited, 0, nAtoms * sizeof(char));
|
|
unsigned count = 1;
|
|
std::vector<int> nextLevelNbrs;
|
|
char *lastLevelNbrs = (char *)malloc(nAtoms * sizeof(char));
|
|
memset(lastLevelNbrs, 0, nAtoms * sizeof(char));
|
|
char *currentLevelNbrs = (char *)malloc(nAtoms * sizeof(char));
|
|
memset(currentLevelNbrs, 0, nAtoms * sizeof(char));
|
|
int *revisitedNeighbors = (int *)malloc(nAtoms * sizeof(int));
|
|
memset(revisitedNeighbors, 0, nAtoms * sizeof(int));
|
|
while (!neighbors.empty()) {
|
|
unsigned int numLevelNbrs = 0;
|
|
nextLevelNbrs.resize(0);
|
|
while (!neighbors.empty()) {
|
|
int nidx = neighbors.front();
|
|
neighbors.pop_front();
|
|
const Canon::canon_atom &atom = atoms[nidx];
|
|
if (!ringInfo->isInitialized() || ringInfo->numAtomRings(atom.atom->getIdx()) < 1) {
|
|
continue;
|
|
}
|
|
lastLevelNbrs[nidx] = 1;
|
|
visited[nidx] = 1;
|
|
for (unsigned int j = 0; j < atom.degree; j++) {
|
|
int iidx = atom.nbrIds[j];
|
|
if (!visited[iidx]) {
|
|
currentLevelNbrs[iidx] = 1;
|
|
numLevelNbrs++;
|
|
visited[iidx] = 1;
|
|
nextLevelNbrs.push_back(iidx);
|
|
}
|
|
}
|
|
}
|
|
for (unsigned i = 0; i < nAtoms; ++i) {
|
|
if (currentLevelNbrs[i]) {
|
|
const Canon::canon_atom &natom = atoms[i];
|
|
for (unsigned int k = 0; k < natom.degree; k++) {
|
|
int jidx = natom.nbrIds[k];
|
|
if (currentLevelNbrs[jidx] || lastLevelNbrs[jidx]) {
|
|
revisitedNeighbors[jidx] += 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
memset(lastLevelNbrs, 0, nAtoms * sizeof(char));
|
|
for (unsigned i = 0; i < nAtoms; ++i) {
|
|
if (currentLevelNbrs[i]) {
|
|
lastLevelNbrs[i] = 1;
|
|
}
|
|
}
|
|
memset(currentLevelNbrs, 0, nAtoms * sizeof(char));
|
|
std::vector<int> tmp;
|
|
tmp.reserve(30);
|
|
for (unsigned i = 0; i < nAtoms; ++i) {
|
|
if (revisitedNeighbors[i] > 0) {
|
|
tmp.push_back(revisitedNeighbors[i]);
|
|
}
|
|
}
|
|
std::sort(tmp.begin(), tmp.end());
|
|
tmp.push_back(-1);
|
|
for (int i : tmp) {
|
|
if (currentRNIdx >= atoms[idx].revistedNeighbors.size()) {
|
|
atoms[idx].revistedNeighbors.resize(
|
|
atoms[idx].revistedNeighbors.size() + 1000);
|
|
}
|
|
atoms[idx].revistedNeighbors[currentRNIdx] = i;
|
|
currentRNIdx++;
|
|
}
|
|
memset(revisitedNeighbors, 0, nAtoms * sizeof(int));
|
|
|
|
atoms[idx].neighborNum.push_back(numLevelNbrs);
|
|
atoms[idx].neighborNum.push_back(-1);
|
|
|
|
neighbors.insert(neighbors.end(), nextLevelNbrs.begin(),
|
|
nextLevelNbrs.end());
|
|
count++;
|
|
}
|
|
atoms[idx].revistedNeighbors.resize(currentRNIdx);
|
|
|
|
free(visited);
|
|
free(currentLevelNbrs);
|
|
free(lastLevelNbrs);
|
|
free(revisitedNeighbors);
|
|
}
|
|
}
|
|
|
|
template <typename T>
|
|
void rankWithFunctor(T &ftor, bool breakTies, int *order,
|
|
bool useSpecial = false, bool useChirality = false,
|
|
const boost::dynamic_bitset<> *atomsInPlay = nullptr,
|
|
const boost::dynamic_bitset<> *bondsInPlay = nullptr) {
|
|
const ROMol &mol = *ftor.dp_mol;
|
|
canon_atom *atoms = ftor.dp_atoms;
|
|
unsigned int nAts = mol.getNumAtoms();
|
|
int *count = (int *)malloc(nAts * sizeof(int));
|
|
int activeset;
|
|
int *next = (int *)malloc(nAts * sizeof(int));
|
|
int *changed = (int *)malloc(nAts * sizeof(int));
|
|
char *touched = (char *)malloc(nAts * sizeof(char));
|
|
memset(touched, 0, nAts * sizeof(char));
|
|
memset(changed, 1, nAts * sizeof(int));
|
|
CreateSinglePartition(nAts, order, count, atoms);
|
|
// ActivatePartitions(nAts,order,count,activeset,next,changed);
|
|
// RefinePartitions(mol,atoms,ftor,false,order,count,activeset,next,changed,touched);
|
|
#ifdef VERBOSE_CANON
|
|
std::cerr << "1--------" << std::endl;
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
std::cerr << order[i] + 1 << " "
|
|
<< " index: " << atoms[order[i]].index
|
|
<< " count: " << count[order[i]] << std::endl;
|
|
}
|
|
#endif
|
|
ftor.df_useNbrs = true;
|
|
ActivatePartitions(nAts, order, count, activeset, next, changed);
|
|
#ifdef VERBOSE_CANON
|
|
std::cerr << "1a--------" << std::endl;
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
std::cerr << order[i] + 1 << " "
|
|
<< " index: " << atoms[order[i]].index
|
|
<< " count: " << count[order[i]] << std::endl;
|
|
}
|
|
#endif
|
|
RefinePartitions(mol, atoms, ftor, true, order, count, activeset, next,
|
|
changed, touched);
|
|
#ifdef VERBOSE_CANON
|
|
std::cerr << "2--------" << std::endl;
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
std::cerr << order[i] + 1 << " "
|
|
<< " index: " << atoms[order[i]].index
|
|
<< " count: " << count[order[i]] << std::endl;
|
|
}
|
|
#endif
|
|
bool ties = false;
|
|
for (unsigned i = 0; i < nAts; ++i) {
|
|
if (!count[i]) {
|
|
ties = true;
|
|
}
|
|
}
|
|
if (useChirality && ties) {
|
|
SpecialChiralityAtomCompareFunctor scftor(atoms, mol, atomsInPlay,
|
|
bondsInPlay);
|
|
ActivatePartitions(nAts, order, count, activeset, next, changed);
|
|
RefinePartitions(mol, atoms, scftor, true, order, count, activeset, next,
|
|
changed, touched);
|
|
#ifdef VERBOSE_CANON
|
|
std::cerr << "2a--------" << std::endl;
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
std::cerr << order[i] + 1 << " "
|
|
<< " index: " << atoms[order[i]].index
|
|
<< " count: " << count[order[i]] << std::endl;
|
|
}
|
|
#endif
|
|
}
|
|
ties = false;
|
|
unsigned symRingAtoms = 0;
|
|
unsigned ringAtoms = 0;
|
|
bool branchingRingAtom = false;
|
|
RingInfo *ringInfo = mol.getRingInfo();
|
|
for (unsigned i = 0; i < nAts; ++i) {
|
|
if (ringInfo->isInitialized() && ringInfo->numAtomRings(order[i])) {
|
|
if (count[order[i]] > 2) {
|
|
symRingAtoms += count[order[i]];
|
|
}
|
|
ringAtoms++;
|
|
if (ringInfo->isInitialized() && ringInfo->numAtomRings(order[i]) > 1 && count[order[i]] > 1) {
|
|
branchingRingAtom = true;
|
|
}
|
|
}
|
|
if (!count[i]) {
|
|
ties = true;
|
|
}
|
|
}
|
|
// std::cout << " " << ringAtoms << " " << symRingAtoms << std::endl;
|
|
if (useSpecial && ties && ringAtoms > 0 &&
|
|
static_cast<float>(symRingAtoms) / ringAtoms > 0.5 && branchingRingAtom) {
|
|
SpecialSymmetryAtomCompareFunctor sftor(atoms, mol, atomsInPlay,
|
|
bondsInPlay);
|
|
compareRingAtomsConcerningNumNeighbors(atoms, nAts, mol);
|
|
ActivatePartitions(nAts, order, count, activeset, next, changed);
|
|
RefinePartitions(mol, atoms, sftor, true, order, count, activeset, next,
|
|
changed, touched);
|
|
#ifdef VERBOSE_CANON
|
|
std::cerr << "2b--------" << std::endl;
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
std::cerr << order[i] + 1 << " "
|
|
<< " index: " << atoms[order[i]].index
|
|
<< " count: " << count[order[i]] << std::endl;
|
|
}
|
|
#endif
|
|
}
|
|
if (breakTies) {
|
|
BreakTies(mol, atoms, ftor, true, order, count, activeset, next, changed,
|
|
touched);
|
|
#ifdef VERBOSE_CANON
|
|
std::cerr << "3--------" << std::endl;
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
std::cerr << order[i] + 1 << " "
|
|
<< " index: " << atoms[order[i]].index
|
|
<< " count: " << count[order[i]] << std::endl;
|
|
}
|
|
#endif
|
|
}
|
|
free(count);
|
|
free(next);
|
|
free(touched);
|
|
free(changed);
|
|
}
|
|
|
|
namespace {
|
|
bool hasRingNbr(const ROMol &mol, const Atom *at) {
|
|
ROMol::ADJ_ITER beg, end;
|
|
boost::tie(beg, end) = mol.getAtomNeighbors(at);
|
|
while (beg != end) {
|
|
const Atom* nbr = mol[*beg];
|
|
++beg;
|
|
if ((nbr->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW ||
|
|
nbr->getChiralTag() == Atom::CHI_TETRAHEDRAL_CCW) &&
|
|
nbr->hasProp(common_properties::_ringStereoAtoms)) {
|
|
return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
void getNbrs(const ROMol &mol, const Atom *at, int *ids) {
|
|
ROMol::OEDGE_ITER beg, end;
|
|
boost::tie(beg, end) = mol.getAtomBonds(at);
|
|
unsigned int idx = 0;
|
|
|
|
while (beg != end) {
|
|
const Bond* bond = (mol)[*beg];
|
|
++beg;
|
|
unsigned int nbrIdx = bond->getOtherAtomIdx(at->getIdx());
|
|
ids[idx] = nbrIdx;
|
|
++idx;
|
|
}
|
|
}
|
|
|
|
bondholder makeBondHolder(const Bond *bond, unsigned int otherIdx,
|
|
bool includeChirality) {
|
|
Bond::BondStereo stereo = Bond::STEREONONE;
|
|
if (includeChirality) {
|
|
stereo = bond->getStereo();
|
|
if (stereo == Bond::STEREOANY) {
|
|
stereo = Bond::STEREONONE;
|
|
}
|
|
}
|
|
Bond::BondType bt =
|
|
bond->getIsAromatic() ? Bond::AROMATIC : bond->getBondType();
|
|
return bondholder(bt, stereo, otherIdx, 0);
|
|
}
|
|
void getBonds(const ROMol &mol, const Atom *at, std::vector<bondholder> &nbrs,
|
|
bool includeChirality) {
|
|
ROMol::OEDGE_ITER beg, end;
|
|
boost::tie(beg, end) = mol.getAtomBonds(at);
|
|
while (beg != end) {
|
|
const Bond* bond = (mol)[*beg];
|
|
++beg;
|
|
nbrs.push_back(makeBondHolder(
|
|
bond, bond->getOtherAtomIdx(at->getIdx()), includeChirality));
|
|
}
|
|
std::sort(nbrs.begin(), nbrs.end(), bondholder::greater);
|
|
}
|
|
|
|
void getChiralBonds(const ROMol &mol, const Atom *at,
|
|
std::vector<bondholder> &nbrs) {
|
|
ROMol::OEDGE_ITER beg, end;
|
|
boost::tie(beg, end) = mol.getAtomBonds(at);
|
|
while (beg != end) {
|
|
const Bond* bond = (mol)[*beg];
|
|
++beg;
|
|
unsigned int nbrIdx = bond->getOtherAtomIdx(at->getIdx());
|
|
const Atom *nbr = mol.getAtomWithIdx(nbrIdx);
|
|
unsigned int degreeNbr = nbr->getDegree();
|
|
unsigned int nReps = 1;
|
|
unsigned int stereo = 0;
|
|
// FIX: Since the user can set the stereoatoms, the use of STEREOCIS and
|
|
// STEREOTRANS here isn't actually canonical. In order for that to work we
|
|
// would need to be able to canonicalize the CIS/TRANS assignment, which is
|
|
// not currently being done
|
|
switch (bond->getStereo()) {
|
|
case Bond::STEREOZ:
|
|
case Bond::STEREOCIS:
|
|
stereo = 1;
|
|
break;
|
|
case Bond::STEREOE:
|
|
case Bond::STEREOTRANS:
|
|
stereo = 2;
|
|
break;
|
|
default:
|
|
stereo = 0;
|
|
}
|
|
if (bond->getBondType() == Bond::DOUBLE && nbr->getAtomicNum() == 15 &&
|
|
(degreeNbr == 4 || degreeNbr == 3)) {
|
|
// a special case for chiral phophorous compounds
|
|
// (this was leading to incorrect assignment of
|
|
// R/S labels ):
|
|
nReps = 1;
|
|
// 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) Contibutions by d orbitals to bonds of quadriligant atoms are
|
|
// neglected."
|
|
// FIX: this applies to more than just P
|
|
} else {
|
|
nReps =
|
|
static_cast<unsigned int>(floor(2. * bond->getBondTypeAsDouble()));
|
|
}
|
|
unsigned int symclass =
|
|
nbr->getAtomicNum() * ATNUM_CLASS_OFFSET + nbrIdx + 1;
|
|
bondholder bh(bondholder(Bond::SINGLE, stereo, nbrIdx, symclass));
|
|
auto iPos = std::lower_bound(nbrs.begin(), nbrs.end(), bh);
|
|
nbrs.insert(iPos, nReps, bh);
|
|
}
|
|
std::reverse(nbrs.begin(), nbrs.end());
|
|
|
|
if (!at->needsUpdatePropertyCache()) {
|
|
for (unsigned int ii = 0; ii < at->getTotalNumHs(); ++ii) {
|
|
nbrs.push_back(bondholder(Bond::SINGLE, Bond::STEREONONE,
|
|
ATNUM_CLASS_OFFSET, ATNUM_CLASS_OFFSET));
|
|
nbrs.push_back(bondholder(Bond::SINGLE, Bond::STEREONONE,
|
|
ATNUM_CLASS_OFFSET, ATNUM_CLASS_OFFSET));
|
|
}
|
|
}
|
|
}
|
|
|
|
void basicInitCanonAtom(const ROMol &mol, Canon::canon_atom &atom,
|
|
const int &idx) {
|
|
atom.atom = mol.getAtomWithIdx(idx);
|
|
atom.index = idx;
|
|
atom.p_symbol = nullptr;
|
|
atom.degree = atom.atom->getDegree();
|
|
atom.nbrIds = (int *)malloc(atom.degree * sizeof(int));
|
|
getNbrs(mol, atom.atom, atom.nbrIds);
|
|
}
|
|
|
|
void advancedInitCanonAtom(const ROMol &mol, Canon::canon_atom &atom,
|
|
const int &idx) {
|
|
RDUNUSED_PARAM(idx);
|
|
atom.totalNumHs = atom.atom->getTotalNumHs();
|
|
atom.isRingStereoAtom =
|
|
(atom.atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW ||
|
|
atom.atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CCW) &&
|
|
atom.atom->hasProp(common_properties::_ringStereoAtoms);
|
|
atom.hasRingNbr = hasRingNbr(mol, atom.atom);
|
|
}
|
|
} // end anonymous namespace
|
|
|
|
void initCanonAtoms(const ROMol &mol, std::vector<Canon::canon_atom> &atoms,
|
|
bool includeChirality) {
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
basicInitCanonAtom(mol, atoms[i], i);
|
|
advancedInitCanonAtom(mol, atoms[i], i);
|
|
atoms[i].bonds.reserve(atoms[i].degree);
|
|
getBonds(mol, atoms[i].atom, atoms[i].bonds, includeChirality);
|
|
}
|
|
}
|
|
|
|
void initFragmentCanonAtoms(const ROMol &mol,
|
|
std::vector<Canon::canon_atom> &atoms,
|
|
bool includeChirality,
|
|
const std::vector<std::string> *atomSymbols,
|
|
const boost::dynamic_bitset<> &atomsInPlay,
|
|
const boost::dynamic_bitset<> &bondsInPlay) {
|
|
// start by initializing the atoms
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
atoms[i].atom = mol.getAtomWithIdx(i);
|
|
atoms[i].index = i;
|
|
// we don't care about overall degree, so we start that at zero, and
|
|
// then count the degree in the fragment itself below.
|
|
atoms[i].degree = 0;
|
|
if (atomsInPlay[i]) {
|
|
atoms[i].nbrIds = (int *)calloc(atoms[i].atom->getDegree(), sizeof(int));
|
|
if (atomSymbols) {
|
|
atoms[i].p_symbol = &(*atomSymbols)[i];
|
|
} else {
|
|
atoms[i].p_symbol = nullptr;
|
|
}
|
|
advancedInitCanonAtom(mol, atoms[i], i);
|
|
atoms[i].bonds.reserve(atoms[i].atom->getDegree());
|
|
}
|
|
}
|
|
|
|
// now deal with the bonds in the fragment.
|
|
// these set the atomic degrees and update the neighbor lists
|
|
for (ROMol::ConstBondIterator bI = mol.beginBonds(); bI != mol.endBonds();
|
|
++bI) {
|
|
if (!bondsInPlay[(*bI)->getIdx()] ||
|
|
!atomsInPlay[(*bI)->getBeginAtomIdx()] ||
|
|
!atomsInPlay[(*bI)->getEndAtomIdx()])
|
|
continue;
|
|
Canon::canon_atom &begAt = atoms[(*bI)->getBeginAtomIdx()];
|
|
Canon::canon_atom &endAt = atoms[(*bI)->getEndAtomIdx()];
|
|
begAt.nbrIds[begAt.degree] = (*bI)->getEndAtomIdx();
|
|
endAt.nbrIds[endAt.degree] = (*bI)->getBeginAtomIdx();
|
|
begAt.degree++;
|
|
endAt.degree++;
|
|
begAt.bonds.push_back(
|
|
makeBondHolder(*bI, (*bI)->getEndAtomIdx(), includeChirality));
|
|
endAt.bonds.push_back(
|
|
makeBondHolder(*bI, (*bI)->getBeginAtomIdx(), includeChirality));
|
|
}
|
|
|
|
// and now we can do the last bit for each atom
|
|
for (size_t i = 0; i < mol.getNumAtoms(); ++i) {
|
|
if (!atomsInPlay[i]) continue;
|
|
// this is the fix for github #1567: we let the atom's degree
|
|
// in the original molecule influence its degree in the fragment
|
|
atoms[i].totalNumHs +=
|
|
(mol.getAtomWithIdx(i)->getDegree() - atoms[i].degree);
|
|
|
|
// and sort our list of neighboring bonds
|
|
std::sort(atoms[i].bonds.begin(), atoms[i].bonds.end(),
|
|
bondholder::greater);
|
|
}
|
|
}
|
|
|
|
void initChiralCanonAtoms(const ROMol &mol,
|
|
std::vector<Canon::canon_atom> &atoms) {
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
basicInitCanonAtom(mol, atoms[i], i);
|
|
getChiralBonds(mol, atoms[i].atom, atoms[i].bonds);
|
|
}
|
|
}
|
|
|
|
void freeCanonAtoms(std::vector<Canon::canon_atom> &atoms) {
|
|
for (auto &atom : atoms) {
|
|
if (atom.nbrIds) {
|
|
free(atom.nbrIds);
|
|
atom.nbrIds = nullptr;
|
|
}
|
|
}
|
|
}
|
|
|
|
void updateAtomNeighborIndex(canon_atom *atoms, std::vector<bondholder> &nbrs) {
|
|
for (auto &nbr : nbrs) {
|
|
unsigned nbrIdx = nbr.nbrIdx;
|
|
unsigned newSymClass = atoms[nbrIdx].index;
|
|
nbr.nbrSymClass = newSymClass;
|
|
}
|
|
std::sort(nbrs.begin(), nbrs.end(), bondholder::greater);
|
|
}
|
|
|
|
void updateAtomNeighborNumSwaps(
|
|
canon_atom *atoms, std::vector<bondholder> &nbrs, unsigned int atomIdx,
|
|
std::vector<std::pair<unsigned int, unsigned int>> &result) {
|
|
for (auto &nbr : nbrs) {
|
|
unsigned nbrIdx = nbr.nbrIdx;
|
|
|
|
if (atoms[nbrIdx].atom->getChiralTag() != 0) {
|
|
std::vector<int> ref, probe;
|
|
for (unsigned i = 0; i < atoms[nbrIdx].degree; ++i) {
|
|
ref.push_back(atoms[nbrIdx].nbrIds[i]);
|
|
}
|
|
|
|
probe.push_back(atomIdx);
|
|
for (auto &bond : atoms[nbrIdx].bonds) {
|
|
if (bond.nbrIdx != atomIdx) {
|
|
probe.push_back(bond.nbrIdx);
|
|
}
|
|
}
|
|
|
|
int nSwaps = static_cast<int>(countSwapsToInterconvert(ref, probe));
|
|
if (atoms[nbrIdx].atom->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW) {
|
|
if (nSwaps % 2) {
|
|
result.push_back(std::make_pair(nbr.nbrSymClass, 2));
|
|
} else {
|
|
result.push_back(std::make_pair(nbr.nbrSymClass, 1));
|
|
}
|
|
} else if (atoms[nbrIdx].atom->getChiralTag() ==
|
|
Atom::CHI_TETRAHEDRAL_CCW) {
|
|
if (nSwaps % 2) {
|
|
result.push_back(std::make_pair(nbr.nbrSymClass, 1));
|
|
} else {
|
|
result.push_back(std::make_pair(nbr.nbrSymClass, 2));
|
|
}
|
|
}
|
|
} else {
|
|
result.push_back(std::make_pair(nbr.nbrSymClass, 0));
|
|
}
|
|
}
|
|
sort(result.begin(), result.end());
|
|
}
|
|
|
|
void rankMolAtoms(const ROMol &mol, std::vector<unsigned int> &res,
|
|
bool breakTies, bool includeChirality, bool includeIsotopes) {
|
|
res.resize(mol.getNumAtoms());
|
|
|
|
if (!mol.getNumAtoms())
|
|
return;
|
|
|
|
std::vector<Canon::canon_atom> atoms(mol.getNumAtoms());
|
|
initCanonAtoms(mol, atoms, includeChirality);
|
|
AtomCompareFunctor ftor(&atoms.front(), mol);
|
|
ftor.df_useIsotopes = includeIsotopes;
|
|
ftor.df_useChirality = includeChirality;
|
|
ftor.df_useChiralityRings = includeChirality;
|
|
|
|
int *order = (int *)malloc(mol.getNumAtoms() * sizeof(int));
|
|
rankWithFunctor(ftor, breakTies, order, true, includeChirality);
|
|
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
res[order[i]] = atoms[order[i]].index;
|
|
}
|
|
free(order);
|
|
freeCanonAtoms(atoms);
|
|
} // end of rankMolAtoms()
|
|
|
|
void rankFragmentAtoms(const ROMol &mol, std::vector<unsigned int> &res,
|
|
const boost::dynamic_bitset<> &atomsInPlay,
|
|
const boost::dynamic_bitset<> &bondsInPlay,
|
|
const std::vector<std::string> *atomSymbols,
|
|
bool breakTies, bool includeChirality,
|
|
bool includeIsotopes) {
|
|
PRECONDITION(atomsInPlay.size() == mol.getNumAtoms(), "bad atomsInPlay size");
|
|
PRECONDITION(bondsInPlay.size() == mol.getNumBonds(), "bad bondsInPlay size");
|
|
PRECONDITION(!atomSymbols || atomSymbols->size() == mol.getNumAtoms(),
|
|
"bad atomSymbols size");
|
|
|
|
res.resize(mol.getNumAtoms());
|
|
|
|
if (!mol.getNumAtoms())
|
|
return;
|
|
|
|
std::vector<Canon::canon_atom> atoms(mol.getNumAtoms());
|
|
initFragmentCanonAtoms(mol, atoms, includeChirality, atomSymbols, atomsInPlay,
|
|
bondsInPlay);
|
|
|
|
AtomCompareFunctor ftor(&atoms.front(), mol, &atomsInPlay, &bondsInPlay);
|
|
ftor.df_useIsotopes = includeIsotopes;
|
|
ftor.df_useChirality = includeChirality;
|
|
|
|
int *order = (int *)malloc(mol.getNumAtoms() * sizeof(int));
|
|
|
|
rankWithFunctor(ftor, breakTies, order, true, includeChirality, &atomsInPlay,
|
|
&bondsInPlay);
|
|
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
res[order[i]] = atoms[order[i]].index;
|
|
}
|
|
free(order);
|
|
freeCanonAtoms(atoms);
|
|
} // end of rankFragmentAtoms()
|
|
|
|
void chiralRankMolAtoms(const ROMol &mol, std::vector<unsigned int> &res) {
|
|
res.resize(mol.getNumAtoms());
|
|
|
|
if(!mol.getNumAtoms())
|
|
return;
|
|
|
|
std::vector<Canon::canon_atom> atoms(mol.getNumAtoms());
|
|
initChiralCanonAtoms(mol, atoms);
|
|
ChiralAtomCompareFunctor ftor(&atoms.front(), mol);
|
|
|
|
int *order = (int *)malloc(mol.getNumAtoms() * sizeof(int));
|
|
rankWithFunctor(ftor, false, order);
|
|
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
res[order[i]] = atoms[order[i]].index;
|
|
}
|
|
free(order);
|
|
freeCanonAtoms(atoms);
|
|
} // end of rankMolAtoms()
|
|
} // end of Canon namespace
|
|
} // end of RDKit namespace
|