From 00a9dc49f7de5e69199526d16ef242e5001e72c4 Mon Sep 17 00:00:00 2001 From: Greg Landrum Date: Fri, 6 Jun 2025 15:30:07 +0200 Subject: [PATCH] Fixes #8559 (#8560) * Fixes #8559 Also includes some minor refactoring of set14Bounds() * apply the same change to the cis/trans trackers * cleanup --------- Co-authored-by: = <=> --- .../DistGeomHelpers/BoundsMatrixBuilder.cpp | 295 +++++++++--------- Code/GraphMol/DistGeomHelpers/catch_tests.cpp | 12 + 2 files changed, 162 insertions(+), 145 deletions(-) diff --git a/Code/GraphMol/DistGeomHelpers/BoundsMatrixBuilder.cpp b/Code/GraphMol/DistGeomHelpers/BoundsMatrixBuilder.cpp index 3612631b1..b683df84b 100644 --- a/Code/GraphMol/DistGeomHelpers/BoundsMatrixBuilder.cpp +++ b/Code/GraphMol/DistGeomHelpers/BoundsMatrixBuilder.cpp @@ -1,5 +1,5 @@ // -// Copyright (C) 2004-2022 Greg Landrum and other RDKit contributors +// Copyright (C) 2004-2025 Greg Landrum and other RDKit contributors // // @@ All Rights Reserved @@ // This file is part of the RDKit. @@ -22,6 +22,7 @@ #include #include #include +#include const double DIST12_DELTA = 0.01; // const double ANGLE_DELTA = 0.0837; @@ -73,11 +74,6 @@ class ComputedData { bondAdj.reset(bAdj); auto *bAngles = new RDNumeric::DoubleSymmMatrix(nBonds, -1.0); bondAngles.reset(bAngles); - DGeomHelpers::BIT_SET::size_type capacity = - boost::numeric_cast( - static_cast(nBonds) * nBonds * nBonds); - cisPaths.resize(capacity); - transPaths.resize(capacity); set15Atoms.resize(nAtoms * nAtoms); } @@ -87,8 +83,8 @@ class ComputedData { SymmIntMatPtr bondAdj; // bond adjacency matrix SymmDoubleMatPtr bondAngles; PATH14_VECT paths14; - BIT_SET cisPaths; - BIT_SET transPaths; + std::unordered_set cisPaths; + std::unordered_set transPaths; BIT_SET set15Atoms; }; @@ -677,18 +673,18 @@ void _setInRing14Bounds(const ROMol &mol, const Bond *bnd1, const Bond *bnd2, dl = RDGeom::compute14DistCis(bl1, bl2, bl3, ba12, ba23) - GEN_DIST_TOL; du = dl + 2 * GEN_DIST_TOL; path14.type = Path14Configuration::CIS; - accumData.cisPaths[static_cast(bid1) * nb * nb + bid2 * nb + - bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + bid2 * nb + - bid1] = 1; + accumData.cisPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.cisPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); } else if (preferTrans) { dl = RDGeom::compute14DistTrans(bl1, bl2, bl3, ba12, ba23) - GEN_DIST_TOL; du = dl + 2 * GEN_DIST_TOL; path14.type = Path14Configuration::TRANS; - accumData.transPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.transPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.transPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.transPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); } else { // basically we will assume 0 to 180 allowed @@ -770,10 +766,10 @@ void _setTwoInSameRing14Bounds(const ROMol &mol, const Bond *bnd1, dl -= GEN_DIST_TOL; du += GEN_DIST_TOL; path14.type = Path14Configuration::TRANS; - accumData.transPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.transPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.transPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.transPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); } else { // here we will assume anything is possible @@ -1015,10 +1011,10 @@ void _setChain14Bounds(const ROMol &mol, const Bond *bnd1, const Bond *bnd2, dl = RDGeom::compute14DistCis(bl1, bl2, bl3, ba12, ba23) - GEN_DIST_TOL; du = dl + 2 * GEN_DIST_TOL; path14.type = Path14Configuration::CIS; - accumData.cisPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.cisPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.cisPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); // BOOST_LOG(rdDebugLog) << "Special 5 " << aid1 << " " << aid4 << // "\n"; } else if (bnd2->getStereo() > Bond::STEREOANY) { @@ -1031,10 +1027,10 @@ void _setChain14Bounds(const ROMol &mol, const Bond *bnd1, const Bond *bnd2, // BOOST_LOG(rdDebugLog) << "Special 6 " << aid1 << " " << aid4 // << // "\n"; - accumData.cisPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.cisPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.cisPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); } else { // BOOST_LOG(rdDebugLog) << "Special 7 " << aid1 << " " << aid4 << // "\n"; @@ -1043,10 +1039,10 @@ void _setChain14Bounds(const ROMol &mol, const Bond *bnd1, const Bond *bnd2, dl -= GEN_DIST_TOL; du += GEN_DIST_TOL; path14.type = Path14Configuration::TRANS; - accumData.transPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.transPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.transPaths.insert( + static_cast(bid1) * nb * nb + bid2 * nb + bid3); + accumData.transPaths.insert( + static_cast(bid3) * nb * nb + bid2 * nb + bid1); } } else { // double bond with no stereo setting can be 0 or 180 @@ -1105,17 +1101,17 @@ void _setChain14Bounds(const ROMol &mol, const Bond *bnd1, const Bond *bnd2, // secondary amide, this is the H, it should be trans to the O dl = RDGeom::compute14DistTrans(bl1, bl2, bl3, ba12, ba23); path14.type = Path14Configuration::TRANS; - accumData.transPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.transPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.transPaths.insert( + static_cast(bid1) * nb * nb + bid2 * nb + bid3); + accumData.transPaths.insert( + static_cast(bid3) * nb * nb + bid2 * nb + bid1); } else { dl = RDGeom::compute14DistCis(bl1, bl2, bl3, ba12, ba23); path14.type = Path14Configuration::CIS; - accumData.cisPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.cisPaths.insert( + static_cast(bid1) * nb * nb + bid2 * nb + bid3); + accumData.cisPaths.insert( + static_cast(bid3) * nb * nb + bid2 * nb + bid1); } du = dl; dl -= GEN_DIST_TOL; @@ -1151,17 +1147,17 @@ void _setChain14Bounds(const ROMol &mol, const Bond *bnd1, const Bond *bnd2, // secondary amide, this is the H, it's cis to atom 5 dl = RDGeom::compute14DistCis(bl1, bl2, bl3, ba12, ba23); path14.type = Path14Configuration::CIS; - accumData.cisPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.cisPaths.insert( + static_cast(bid1) * nb * nb + bid2 * nb + bid3); + accumData.cisPaths.insert( + static_cast(bid3) * nb * nb + bid2 * nb + bid1); } else { dl = RDGeom::compute14DistTrans(bl1, bl2, bl3, ba12, ba23); path14.type = Path14Configuration::TRANS; - accumData.transPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.transPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.transPaths.insert( + static_cast(bid1) * nb * nb + bid2 * nb + bid3); + accumData.transPaths.insert( + static_cast(bid3) * nb * nb + bid2 * nb + bid1); } du = dl; dl -= GEN_DIST_TOL; @@ -1208,10 +1204,10 @@ void _record14Path(const ROMol &mol, unsigned int bid1, unsigned int bid2, path14.bid3 = bid3; if ((ahyb2 == Atom::SP2) && (ahyb3 == Atom::SP2)) { // FIX: check for trans path14.type = Path14Configuration::CIS; - accumData.cisPaths[static_cast(bid1) * nb * nb + bid2 * nb + - bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + bid2 * nb + - bid1] = 1; + accumData.cisPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.cisPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); } else { path14.type = Path14Configuration::OTHER; } @@ -1305,8 +1301,10 @@ void _setMacrocycleTwoInSameRing14Bounds(const ROMol &mol, const Bond *bnd1, atm1))) { dl = RDGeom::compute14DistCis(bl1, bl2, bl3, ba12, ba23); path14.type = Path14Configuration::CIS; - accumData.cisPaths[bid1 * nb * nb + bid2 * nb + bid3] = 1; - accumData.cisPaths[bid3 * nb * nb + bid2 * nb + bid1] = 1; + accumData.cisPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.cisPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); du = dl; dl -= GEN_DIST_TOL; du += GEN_DIST_TOL; @@ -1384,10 +1382,10 @@ void _setMacrocycleAllInSameRing14Bounds(const ROMol &mol, const Bond *bnd1, dl = RDGeom::compute14DistCis(bl1, bl2, bl3, ba12, ba23) - GEN_DIST_TOL; du = dl + 2 * GEN_DIST_TOL; path14.type = Path14Configuration::CIS; - accumData.cisPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.cisPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.cisPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); // BOOST_LOG(rdDebugLog) << "Special 5 " << aid1 << " " << aid4 << // "\n"; } else if (bnd2->getStereo() > Bond::STEREOANY) { @@ -1399,10 +1397,10 @@ void _setMacrocycleAllInSameRing14Bounds(const ROMol &mol, const Bond *bnd1, path14.type = Path14Configuration::CIS; // BOOST_LOG(rdDebugLog) << "Special 6 " << aid1 << " " << aid4 << // "\n"; - accumData.cisPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.cisPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.cisPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); } else { // BOOST_LOG(rdDebugLog) << "Special 7 " << aid1 << " " << aid4 << // "\n"; @@ -1411,10 +1409,10 @@ void _setMacrocycleAllInSameRing14Bounds(const ROMol &mol, const Bond *bnd1, dl -= GEN_DIST_TOL; du += GEN_DIST_TOL; path14.type = Path14Configuration::TRANS; - accumData.transPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.transPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.transPaths.insert( + static_cast(bid1) * nb * nb + bid2 * nb + bid3); + accumData.transPaths.insert( + static_cast(bid3) * nb * nb + bid2 * nb + bid1); } } else { // double bond with no stereo setting can be 0 or 180 @@ -1448,10 +1446,10 @@ void _setMacrocycleAllInSameRing14Bounds(const ROMol &mol, const Bond *bnd1, // still a bit too short, thus we add an additional 0.1, which // is the max that works without triangular smoothing error path14.type = Path14Configuration::TRANS; - accumData.transPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.transPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.transPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.transPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); du = dl; dl -= GEN_DIST_TOL; @@ -1467,10 +1465,10 @@ void _setMacrocycleAllInSameRing14Bounds(const ROMol &mol, const Bond *bnd1, // amide is trans, we're cis: dl = RDGeom::compute14DistCis(bl1, bl2, bl3, ba12, ba23); path14.type = Path14Configuration::CIS; - accumData.cisPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.cisPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.cisPaths.insert(static_cast(bid1) * nb * nb + + bid2 * nb + bid3); + accumData.cisPaths.insert(static_cast(bid3) * nb * nb + + bid2 * nb + bid1); #else // amide is cis, we're trans: if (atm2->getAtomicNum() == 7 && atm2->getDegree() == 3 && @@ -1480,10 +1478,10 @@ void _setMacrocycleAllInSameRing14Bounds(const ROMol &mol, const Bond *bnd1, } else { dl = RDGeom::compute14DistTrans(bl1, bl2, bl3, ba12, ba23); path14.type = Path14Configuration::TRANS; - accumData.transPaths[static_cast(bid1) * nb * nb + - bid2 * nb + bid3] = 1; - accumData.transPaths[static_cast(bid3) * nb * nb + - bid2 * nb + bid1] = 1; + accumData.transPaths.insert( + static_cast(bid1) * nb * nb + bid2 * nb + bid3); + accumData.transPaths.insert( + static_cast(bid3) * nb * nb + bid2 * nb + bid1); } #endif du = dl; @@ -1521,43 +1519,41 @@ void set14Bounds(const ROMol &mol, DistGeom::BoundsMatPtr mmat, bool useMacrocycle14config, bool forceTransAmides) { unsigned int npt = mmat->numRows(); CHECK_INVARIANT(npt == mol.getNumAtoms(), "Wrong size metric matrix"); - - const RingInfo *rinfo = - mol.getRingInfo(); // FIX: make sure we have ring info + // this is 2.6 million bonds, so it's extremly unlikely to ever occur, but we + // might as well check: + const size_t MAX_NUM_BONDS = static_cast( + std::pow(std::numeric_limits::max(), 1. / 3)); + if (mol.getNumBonds() >= MAX_NUM_BONDS) { + throw ValueErrorException( + "Too many bonds in the molecule, cannot compute 1-4 bounds"); + } + const auto rinfo = mol.getRingInfo(); // FIX: make sure we have ring info CHECK_INVARIANT(rinfo, ""); - const VECT_INT_VECT &bondRings = rinfo->bondRings(); - VECT_INT_VECT_CI rii; - unsigned int i, aid2, aid3; - unsigned int bid1, bid2, bid3; - ROMol::OEDGE_ITER beg1, beg2, end1, end2; - ROMol::ConstBondIterator bi; + const auto &bondRings = rinfo->bondRings(); - std::set bidIsMacrocycle; - unsigned int nb = mol.getNumBonds(); - BIT_SET ringBondPairs(nb * nb), donePaths(nb * nb * nb); - unsigned int id1, id2, pid1, pid2, pid3, pid4; + std::unordered_set bidIsMacrocycle; + std::unordered_set ringBondPairs; + std::unordered_set donePaths; + std::uint64_t nb = mol.getNumBonds(); // first we will deal with 1-4 atoms that belong to the same ring - for (rii = bondRings.begin(); rii != bondRings.end(); rii++) { - // we don't need deal with 3 membered rings - unsigned int rSize = rii->size(); + for (const auto &bring : bondRings) { + const auto rSize = bring.size(); + if (rSize < 3) { + continue; // rings with less than 3 bonds are not useful + } + auto bid1 = bring[rSize - 1]; + for (auto i = 0u; i < rSize; i++) { + auto bid2 = bring[i]; + auto bid3 = bring[(i + 1) % rSize]; + auto pid1 = bid1 * nb + bid2; + auto pid2 = bid2 * nb + bid1; + auto id1 = bid1 * nb * nb + bid2 * nb + bid3; + auto id2 = bid3 * nb * nb + bid2 * nb + bid1; - bid1 = (*rii)[rSize - 1]; - for (i = 0; i < rSize; i++) { - bid2 = (*rii)[i]; - if (i == rSize - 1) { - bid3 = (*rii)[0]; - } else { - bid3 = (*rii)[i + 1]; - } - pid1 = bid1 * nb + bid2; - pid2 = bid2 * nb + bid1; - id1 = bid1 * nb * nb + bid2 * nb + bid3; - id2 = bid3 * nb * nb + bid2 * nb + bid1; - - ringBondPairs[pid1] = 1; - ringBondPairs[pid2] = 1; - donePaths[id1] = 1; - donePaths[id2] = 1; + ringBondPairs.insert(pid1); + ringBondPairs.insert(pid2); + donePaths.insert(id1); + donePaths.insert(id2); if (rSize > 5) { if (useMacrocycle14config && rSize >= minMacrocycleRingSize) { @@ -1577,40 +1573,39 @@ void set14Bounds(const ROMol &mol, DistGeom::BoundsMatPtr mmat, bid1 = bid2; } // loop over bonds in the ring } // end of all rings - for (bi = mol.beginBonds(); bi != mol.endBonds(); bi++) { - bid2 = (*bi)->getIdx(); - aid2 = (*bi)->getBeginAtomIdx(); - aid3 = (*bi)->getEndAtomIdx(); - boost::tie(beg1, end1) = mol.getAtomBonds(mol.getAtomWithIdx(aid2)); - while (beg1 != end1) { - const Bond *bnd1 = mol[*beg1]; - bid1 = bnd1->getIdx(); + for (const auto bond : mol.bonds()) { + auto bid2 = bond->getIdx(); + auto aid2 = bond->getBeginAtomIdx(); + auto aid3 = bond->getEndAtomIdx(); + for (const auto bnd1 : mol.atomBonds(mol.getAtomWithIdx(aid2))) { + auto bid1 = bnd1->getIdx(); if (bid1 != bid2) { - boost::tie(beg2, end2) = mol.getAtomBonds(mol.getAtomWithIdx(aid3)); - while (beg2 != end2) { - const Bond *bnd3 = mol[*beg2]; - bid3 = bnd3->getIdx(); + for (const auto bnd3 : mol.atomBonds(mol.getAtomWithIdx(aid3))) { + auto bid3 = bnd3->getIdx(); if (bid3 != bid2) { - id1 = nb * nb * bid1 + nb * bid2 + bid3; - id2 = nb * nb * bid3 + nb * bid2 + bid1; - if ((!donePaths[id1]) && (!donePaths[id2])) { + auto id1 = bid1 * nb * nb + bid2 * nb + bid3; + auto id2 = bid3 * nb * nb + bid2 * nb + bid1; + if (donePaths.find(id1) == donePaths.end() && + donePaths.find(id2) == donePaths.end()) { // we haven't dealt with this path before - pid1 = bid1 * nb + bid2; - pid2 = bid2 * nb + bid1; - pid3 = bid2 * nb + bid3; - pid4 = bid3 * nb + bid2; + auto pid1 = bid1 * nb + bid2; + auto pid2 = bid2 * nb + bid1; + auto pid3 = bid2 * nb + bid3; + auto pid4 = bid3 * nb + bid2; - if (ringBondPairs[pid1] || ringBondPairs[pid2] || - ringBondPairs[pid3] || ringBondPairs[pid4]) { + if (ringBondPairs.find(pid1) != ringBondPairs.end() || + ringBondPairs.find(pid2) != ringBondPairs.end() || + ringBondPairs.find(pid3) != ringBondPairs.end() || + ringBondPairs.find(pid4) != ringBondPairs.end()) { // either (bid1, bid2) or (bid2, bid3) are in the // same ring (note all three cannot be in the same // ring; we dealt with that before) if (useMacrocycle14config && bidIsMacrocycle.find(bid2) != bidIsMacrocycle.end()) { _setMacrocycleTwoInSameRing14Bounds( - mol, bnd1, (*bi), bnd3, accumData, mmat, distMatrix); + mol, bnd1, bond, bnd3, accumData, mmat, distMatrix); } else { - _setTwoInSameRing14Bounds(mol, bnd1, (*bi), bnd3, accumData, + _setTwoInSameRing14Bounds(mol, bnd1, bond, bnd3, accumData, mmat, distMatrix); } } else if (((rinfo->numBondRings(bid1) > 0) && @@ -1625,26 +1620,24 @@ void set14Bounds(const ROMol &mol, DistGeom::BoundsMatPtr mmat, // bid2 are ring bonds that belong to ring r1 and // r2, then bid3 is either an external bond or // belongs to a third ring r3. - _setTwoInDiffRing14Bounds(mol, bnd1, (*bi), bnd3, accumData, + _setTwoInDiffRing14Bounds(mol, bnd1, bond, bnd3, accumData, mmat, distMatrix); } else if (rinfo->numBondRings(bid2) > 0) { // the middle bond is a ring bond and the other // two do not belong to the same ring or are // non-ring bonds - _setShareRingBond14Bounds(mol, bnd1, (*bi), bnd3, accumData, + _setShareRingBond14Bounds(mol, bnd1, bond, bnd3, accumData, mmat, distMatrix); } else { // middle bond not a ring - _setChain14Bounds(mol, bnd1, (*bi), bnd3, accumData, mmat, + _setChain14Bounds(mol, bnd1, bond, bnd3, accumData, mmat, distMatrix, forceTransAmides); } } } - ++beg2; } } - ++beg1; } } } @@ -1674,6 +1667,15 @@ void setTopolBounds(const ROMol &mol, DistGeom::BoundsMatPtr mmat, if (!na) { throw ValueErrorException("molecule has no atoms"); } + // this is 2.6 million bonds, so it's extremly unlikely to ever occur, but we + // might as well check: + const size_t MAX_NUM_BONDS = static_cast( + std::pow(std::numeric_limits::max(), 1. / 3)); + if (mol.getNumBonds() >= MAX_NUM_BONDS) { + throw ValueErrorException( + "Too many bonds in the molecule, cannot compute 1-4 bounds"); + } + ComputedData accumData(na, nb); double *distMatrix = nullptr; distMatrix = MolOps::getDistanceMat(mol); @@ -1990,11 +1992,12 @@ void _set15BoundsHelper(const ROMol &mol, unsigned int bid1, unsigned int bid2, unsigned long pathId = static_cast(bid2) * nb * nb + (bid3)*nb + i; if (type == 0) { - if (accumData.cisPaths[pathId]) { + if (accumData.cisPaths.find(pathId) != accumData.cisPaths.end()) { dl = _compute15DistsCisCis(d1, d2, d3, d4, ang12, ang23, ang34); du = dl + DIST15_TOL; dl -= DIST15_TOL; - } else if (accumData.transPaths[pathId]) { + } else if (accumData.transPaths.find(pathId) != + accumData.transPaths.end()) { dl = _compute15DistsCisTrans(d1, d2, d3, d4, ang12, ang23, ang34); du = dl + DIST15_TOL; dl -= DIST15_TOL; @@ -2007,11 +2010,12 @@ void _set15BoundsHelper(const ROMol &mol, unsigned int bid1, unsigned int bid2, } } else if (type == 1) { - if (accumData.cisPaths[pathId]) { + if (accumData.cisPaths.find(pathId) != accumData.cisPaths.end()) { dl = _compute15DistsTransCis(d1, d2, d3, d4, ang12, ang23, ang34); du = dl + DIST15_TOL; dl -= DIST15_TOL; - } else if (accumData.transPaths[pathId]) { + } else if (accumData.transPaths.find(pathId) != + accumData.transPaths.end()) { dl = _compute15DistsTransTrans(d1, d2, d3, d4, ang12, ang23, ang34); du = dl + DIST15_TOL; @@ -2025,13 +2029,14 @@ void _set15BoundsHelper(const ROMol &mol, unsigned int bid1, unsigned int bid2, DIST15_TOL; } } else { - if (accumData.cisPaths[pathId]) { + if (accumData.cisPaths.find(pathId) != accumData.cisPaths.end()) { dl = _compute15DistsCisCis(d4, d3, d2, d1, ang34, ang23, ang12) - DIST15_TOL; du = _compute15DistsCisTrans(d4, d3, d2, d1, ang34, ang23, ang12) + DIST15_TOL; - } else if (accumData.transPaths[pathId]) { + } else if (accumData.transPaths.find(pathId) != + accumData.transPaths.end()) { dl = _compute15DistsTransCis(d4, d3, d2, d1, ang34, ang23, ang12) - DIST15_TOL; diff --git a/Code/GraphMol/DistGeomHelpers/catch_tests.cpp b/Code/GraphMol/DistGeomHelpers/catch_tests.cpp index 2239af07a..3252bf581 100644 --- a/Code/GraphMol/DistGeomHelpers/catch_tests.cpp +++ b/Code/GraphMol/DistGeomHelpers/catch_tests.cpp @@ -1234,3 +1234,15 @@ TEST_CASE("FindDoubleBonds") { CHECK(stereoDoubleBonds.empty()); } } + +TEST_CASE("Github #8559: seg fault in setTopolBounds") { + SECTION("as reported") { + auto mol = + "CC1CN([C@H]2C[C@H](O[PH](O)(O)OC[C@H]3O[C@@H](N4CNC5C(N)NCNC54)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CC(C)C(O)NC4O)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CNC5C(N)NCNC54)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CCC(N)NC4O)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CNC5C(O)NC(N)NC54)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CNC5C(N)NCNC54)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CNC5C(N)NCNC54)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CNC5C(O)NC(N)NC54)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CC(C)C(O)NC4O)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CC(C)C(O)NC4O)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CNC5C(N)NCNC54)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CC(C)C(O)NC4O)C[C@@H]3O[PH](O)(O)OC[C@H]3O[C@@H](N4CCC(N)NC4O)C[C@@H]3O)[C@@H](CO[PH](O)(O)O[C@H]3C[C@H](N4CCC(N)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(O)NC(N)NC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CC(C)C(O)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(N)NCNC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CC(C)C(O)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(O)NC(N)NC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CC(C)C(O)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(N)NCNC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(N)NCNC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CC(C)C(O)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(N)NCNC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CC(C)C(O)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(O)NC(N)NC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CCC(N)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CC(C)C(O)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CCC(O)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CCC(N)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(N)NCNC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(N)NCNC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CC(C)C(O)NC4O)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(N)NCNC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CNC5C(O)NC(N)NC54)O[C@@H]3CO[PH](O)(O)O[C@H]3C[C@H](N4CCC(N)NC4O)O[C@@H]3CO)O2)C(O)NC1O"_smiles; + REQUIRE(mol); + MolOps::addHs(*mol); + DistGeom::BoundsMatPtr bm{new DistGeom::BoundsMatrix(mol->getNumAtoms())}; + DGeomHelpers::initBoundsMat(bm, 0.0, 1000.0); + DGeomHelpers::setTopolBounds(*mol, bm); + } +} \ No newline at end of file