diff --git a/Code/ForceField/UFF/Utils.cpp b/Code/ForceField/UFF/Utils.cpp index aff454386..f9cf38528 100644 --- a/Code/ForceField/UFF/Utils.cpp +++ b/Code/ForceField/UFF/Utils.cpp @@ -19,17 +19,23 @@ double calculateCosY(const RDGeom::Point3D &iPoint, const RDGeom::Point3D &jPoint, const RDGeom::Point3D &kPoint, const RDGeom::Point3D &lPoint) { + constexpr double zeroTol = 1.0e-16; RDGeom::Point3D rJI = iPoint - jPoint; RDGeom::Point3D rJK = kPoint - jPoint; RDGeom::Point3D rJL = lPoint - jPoint; - rJI.normalize(); - rJK.normalize(); - rJL.normalize(); - + auto l2JI = rJI.lengthSq(); + auto l2JK = rJK.lengthSq(); + auto l2JL = rJL.lengthSq(); + if (l2JI < zeroTol || l2JK < zeroTol || l2JL < zeroTol) { + return 0.0; + } RDGeom::Point3D n = rJI.crossProduct(rJK); - n.normalize(); - - return n.dotProduct(rJL); + n /= (sqrt(l2JI) * sqrt(l2JK)); + auto l2n = n.lengthSq(); + if (l2n < zeroTol) { + return 0.0; + } + return n.dotProduct(rJL) / (sqrt(l2JL) * sqrt(l2n)); } std::tuple diff --git a/Code/Geometry/point.h b/Code/Geometry/point.h index 0a0d3d5a2..8c91c0925 100644 --- a/Code/Geometry/point.h +++ b/Code/Geometry/point.h @@ -145,6 +145,10 @@ class RDKIT_RDGEOMETRYLIB_EXPORT Point3D : public Point { void normalize() override { double l = this->length(); + if (l < zero_tolerance) { + throw std::runtime_error("Cannot normalize a zero length vector"); + } + x /= l; y /= l; z /= l; @@ -196,7 +200,7 @@ class RDKIT_RDGEOMETRYLIB_EXPORT Point3D : public Point { double signedAngleTo(const Point3D &other) const { double res = this->angleTo(other); // check the sign of the z component of the cross product: - if ((this->x * other.y - this->y * other.x) < -1e-6) { + if ((this->x * other.y - this->y * other.x) < -zero_tolerance) { res = 2.0 * M_PI - res; } return res; @@ -357,6 +361,10 @@ class RDKIT_RDGEOMETRYLIB_EXPORT Point2D : public Point { void normalize() override { double ln = this->length(); + if (ln < zero_tolerance) { + throw std::runtime_error("Cannot normalize a zero length vector"); + } + x /= ln; y /= ln; } @@ -400,7 +408,7 @@ class RDKIT_RDGEOMETRYLIB_EXPORT Point2D : public Point { double signedAngleTo(const Point2D &other) const { double res = this->angleTo(other); - if ((this->x * other.y - this->y * other.x) < -1e-6) { + if ((this->x * other.y - this->y * other.x) < -zero_tolerance) { res = 2.0 * M_PI - res; } return res; diff --git a/Code/GraphMol/AddHs.cpp b/Code/GraphMol/AddHs.cpp index 25b293d46..114e67e58 100644 --- a/Code/GraphMol/AddHs.cpp +++ b/Code/GraphMol/AddHs.cpp @@ -19,6 +19,8 @@ #include #include +constexpr double sq_dist_zero_tol = 1.e-4; + namespace RDKit { // Local utility functionality: @@ -161,7 +163,7 @@ RDGeom::Point3D pickBisector(const RDGeom::Point3D &nbr1Vect, const RDGeom::Point3D &nbr2Vect, const RDGeom::Point3D &nbr3Vect) { auto dirVect = nbr2Vect + nbr3Vect; - if (dirVect.lengthSq() < 1e-4) { + if (dirVect.lengthSq() < sq_dist_zero_tol) { // nbr2Vect and nbr3Vect are anti-parallel (was #3854) dirVect = nbr2Vect; std::swap(dirVect.x, dirVect.y); @@ -230,7 +232,7 @@ void setTerminalAtomCoords(ROMol &mol, unsigned int idx, RDGeom::Point3D nbr1Pos = (*cfi)->getAtomPos(nbr1->getIdx()); // get a normalized vector pointing away from the neighbor: nbr1Vect = nbr1Pos - otherPos; - if (nbr1Vect.lengthSq() < 1e-4) { + if (nbr1Vect.lengthSq() < sq_dist_zero_tol) { // no difference, which likely indicates that we have redundant atoms. // just put it on top of the heavy atom. This was #678 (*cfi)->setAtomPos(idx, otherPos); @@ -320,7 +322,8 @@ void setTerminalAtomCoords(ROMol &mol, unsigned int idx, otherPos = (*cfi)->getAtomPos(otherIdx); nbr1Vect = otherPos - (*cfi)->getAtomPos(nbr1->getIdx()); nbr2Vect = otherPos - (*cfi)->getAtomPos(nbr2->getIdx()); - if (nbr1Vect.lengthSq() < 1e-4 || nbr2Vect.lengthSq() < 1e-4) { + if (nbr1Vect.lengthSq() < sq_dist_zero_tol || + nbr2Vect.lengthSq() < sq_dist_zero_tol) { // no difference, which likely indicates that we have redundant atoms. // just put it on top of the heavy atom. This was #678 (*cfi)->setAtomPos(idx, otherPos); @@ -330,6 +333,11 @@ void setTerminalAtomCoords(ROMol &mol, unsigned int idx, nbr2Vect.normalize(); dirVect = nbr1Vect + nbr2Vect; + if (dirVect.lengthSq() < sq_dist_zero_tol) { + // nbr1Vect and nbr2Vect are non-null, but they may + // still cancel each other out + continue; + } dirVect.normalize(); if ((*cfi)->is3D()) { switch (otherAtom->getHybridization()) { @@ -400,8 +408,9 @@ void setTerminalAtomCoords(ROMol &mol, unsigned int idx, nbr1Vect = otherPos - (*cfi)->getAtomPos(nbr1->getIdx()); nbr2Vect = otherPos - (*cfi)->getAtomPos(nbr2->getIdx()); nbr3Vect = otherPos - (*cfi)->getAtomPos(nbr3->getIdx()); - if (nbr1Vect.lengthSq() < 1e-4 || nbr2Vect.lengthSq() < 1e-4 || - nbr3Vect.lengthSq() < 1e-4) { + if (nbr1Vect.lengthSq() < sq_dist_zero_tol || + nbr2Vect.lengthSq() < sq_dist_zero_tol || + nbr3Vect.lengthSq() < sq_dist_zero_tol) { // no difference, which likely indicates that we have redundant atoms. // just put it on top of the heavy atom. This was #678 (*cfi)->setAtomPos(idx, otherPos); @@ -420,6 +429,23 @@ void setTerminalAtomCoords(ROMol &mol, unsigned int idx, 0.1) { // compute the normal: dirVect = nbr1Vect.crossProduct(nbr2Vect); + + // Each of the nbr vectors is non-null, but there might be pairs + // that cancel each other out. Try to find a direction from atoms + // that do not overlap. + if (dirVect.lengthSq() < sq_dist_zero_tol) { + // This definition of dirVect reverses the parity around otherIdx + // the change of sign restores it + dirVect = nbr1Vect.crossProduct(nbr3Vect) * -1; + } + if (dirVect.lengthSq() < sq_dist_zero_tol) { + dirVect = nbr2Vect.crossProduct(nbr3Vect); + } + // We couldn't find a good direction + if (dirVect.lengthSq() < sq_dist_zero_tol) { + continue; + } + std::string cipCode; if (otherAtom->getPropIfPresent(common_properties::_CIPCode, cipCode)) { diff --git a/Code/GraphMol/CMakeLists.txt b/Code/GraphMol/CMakeLists.txt index 8e1ed337c..3e8c5c6f5 100644 --- a/Code/GraphMol/CMakeLists.txt +++ b/Code/GraphMol/CMakeLists.txt @@ -200,5 +200,5 @@ rdkit_catch_test(tableTestsCatch catch_periodictable.cpp LINK_LIBRARIES GraphMol) rdkit_catch_test(molopsTestsCatch catch_molops.cpp - LINK_LIBRARIES SmilesParse GraphMol) + LINK_LIBRARIES FileParsers SmilesParse GraphMol) diff --git a/Code/GraphMol/WedgeBonds.cpp b/Code/GraphMol/WedgeBonds.cpp index 90062820d..19b064c47 100644 --- a/Code/GraphMol/WedgeBonds.cpp +++ b/Code/GraphMol/WedgeBonds.cpp @@ -145,7 +145,15 @@ Bond::BondDir determineBondWedgeState(const Bond *bond, auto tmpPt = conf->getAtomPos(bondAtom->getIdx()); centerLoc.z = 0.0; tmpPt.z = 0.0; - RDGeom::Point3D refVect = centerLoc.directionVector(tmpPt); + + RDGeom::Point3D refVect; + try { + refVect = centerLoc.directionVector(tmpPt); + } catch (const std::runtime_error &) { + // we have a problem with the reference bond; + // it's probably that the center and the tmp atom overlap + return res; + } neighborBondIndices.push_back(bond->getIdx()); neighborBondAngles.push_back(0.0); @@ -154,7 +162,14 @@ Bond::BondDir determineBondWedgeState(const Bond *bond, if (nbrBond != bond) { tmpPt = conf->getAtomPos(otherAtom->getIdx()); tmpPt.z = 0.0; - auto tmpVect = centerLoc.directionVector(tmpPt); + RDGeom::Point3D tmpVect; + try { + tmpVect = centerLoc.directionVector(tmpPt); + } catch (const std::runtime_error &) { + // we have a problem with the tmp bond; + // it's probably that the atoms overlap + return res; + } auto angle = refVect.signedAngleTo(tmpVect); if (angle < 0.0) { angle += 2. * M_PI; @@ -344,14 +359,12 @@ int pickBondToWedge( std::map> pickBondsToWedge( const ROMol &mol, const BondWedgingParameters *params) { - const Conformer *conf = nullptr; if (mol.getNumConformers()) { conf = &mol.getConformer(); } return pickBondsToWedge(mol, params, conf); - } std::map> pickBondsToWedge( diff --git a/Code/GraphMol/catch_molops.cpp b/Code/GraphMol/catch_molops.cpp index 6db14e4f4..f549d4d6d 100644 --- a/Code/GraphMol/catch_molops.cpp +++ b/Code/GraphMol/catch_molops.cpp @@ -14,6 +14,7 @@ #include #include +#include #include #include @@ -211,3 +212,33 @@ TEST_CASE("getMolFormula") { CHECK(formula == "C2[13C]2H10[2H]BO2"); } } + +TEST_CASE("check division by zero in setTerminalAtomCoords") { + SECTION("degree 4") { + auto m = R"CTAB( + Mrv2311 11162401483D + + 0 0 0 0 0 999 V3000 +M V30 BEGIN CTAB +M V30 COUNTS 5 4 0 0 0 +M V30 BEGIN ATOM +M V30 1 C 0 0 0 0 +M V30 2 F -0.7971 -0.9945 1.275 0 +M V30 3 F 0.7971 0.9945 -1.275 0 +M V30 4 F -1.2069 -0.1568 -0.8768 0 +M V30 5 Cl 1.1223 1.1513 0.8154 0 +M V30 END ATOM +M V30 BEGIN BOND +M V30 1 1 1 2 +M V30 2 1 1 3 +M V30 3 1 1 4 +M V30 4 1 1 5 +M V30 END BOND +M V30 END CTAB +M END +)CTAB"_ctab; + REQUIRE(m); + + CHECK_NOTHROW(MolOps::setTerminalAtomCoords(*m, 4, 0)); + } +} \ No newline at end of file diff --git a/Code/Numerics/Vector.h b/Code/Numerics/Vector.h index 43424ea38..7fb26e6c1 100644 --- a/Code/Numerics/Vector.h +++ b/Code/Numerics/Vector.h @@ -22,6 +22,8 @@ #include #include +static constexpr double zero_tolerance = 1.e-16; + namespace RDNumeric { //! A class to represent vectors of numbers. @@ -260,6 +262,9 @@ class Vector { //! Normalize the vector using the L2 norm inline void normalize() { TYPE val = this->normL2(); + if (val < zero_tolerance) { + throw std::runtime_error("Cannot normalize a zero length vector"); + } (*this) /= val; }