diff --git a/Code/GraphMol/GaussianShape/SingleConformerAlignment.cpp b/Code/GraphMol/GaussianShape/SingleConformerAlignment.cpp index 1f23360a9..edc2bf4e5 100644 --- a/Code/GraphMol/GaussianShape/SingleConformerAlignment.cpp +++ b/Code/GraphMol/GaussianShape/SingleConformerAlignment.cpp @@ -82,13 +82,11 @@ void SingleConformerAlignment::getFinalQuatTrans( -d_initQuatTrans[4], -d_initQuatTrans[5], -d_initQuatTrans[6]}); RDGeom::Transform3D initialRot; initialRot.SetRotationFromQuaternion(d_initQuatTrans.data()); - auto tt = reverseInitialTrans * tmp * initialRot; - xform = tt; + xform = reverseInitialTrans * tmp * initialRot; } -std::array SingleConformerAlignment::calcScores(const double *ref, - const double *fit, - bool includeColor) { +std::array SingleConformerAlignment::calcScores( + const double *ref, const double *fit, const bool includeColor) { std::array scores{0.0, 0.0, 0.0, 0.0, 0.0}; scores[3] = calcVolAndGrads(ref, d_nRefShape, d_refCarbonRadii, fit, d_nFitShape, d_fitCarbonRadii, d_gradConverters, @@ -105,13 +103,15 @@ std::array SingleConformerAlignment::calcScores(const double *ref, return scores; } -std::array SingleConformerAlignment::calcScores(bool includeColor) { +std::array SingleConformerAlignment::calcScores( + const bool includeColor) { applyQuatTrans(d_quatTrans); return calcScores(d_refTemp.data(), d_fitTemp.data(), includeColor); } std::array SingleConformerAlignment::calcScores( - const double shapeOvVol, const double colorOvVol, bool includeColor) const { + const double shapeOvVol, const double colorOvVol, + const bool includeColor) const { std::array scores{0.0, 0.0, 0.0, 0.0, 0.0}; scores[3] = shapeOvVol; scores[4] = colorOvVol; @@ -135,31 +135,32 @@ namespace { // This uses the chain rule: the dV/qQ = (dV/dr) * (dr/dQ) where V is the // volume overlap and r is the Cartesian space. Assumes gradConverters // is already the correct size. -void cartToQuatGrads(const double *quat, const double *mol, int numBPts, - std::vector &gradConverters, int gradConvOffset) { +void cartToQuatGrads(const double *quat, const double *mol, const int numBPts, + std::vector &gradConverters, + const int gradConvOffset) { // for ease of ref - auto q = quat[0]; - auto r = quat[1]; - auto s = quat[2]; - auto u = quat[3]; - auto coef = 1.0 / (q * q + r * r + s * s + u * u); + const auto q = quat[0]; + const auto r = quat[1]; + const auto s = quat[2]; + const auto u = quat[3]; + const auto coef = 1.0 / (q * q + r * r + s * s + u * u); for (int i = 0, j = gradConvOffset, k = 12 * gradConvOffset; i < 4 * numBPts; i += 4, ++j, k += 12) { - auto x = mol[i]; - auto y = mol[i + 1]; - auto z = mol[i + 2]; - auto dx_dq = coef * 2.0 * (q * x + u * y - s * z); - auto dx_dr = coef * 2.0 * (r * x + s * y + u * z); - auto dy_dr = coef * 2.0 * (s * x - r * y + q * z); - auto dx_du = coef * 2.0 * (-u * x + q * y + r * z); - auto dz_ds = dx_dq; - auto dy_du = -dx_dq; - auto dy_ds = dx_dr; - auto dz_du = dx_dr; - auto dx_ds = -dy_dr; - auto dz_dq = dy_dr; - auto dy_dq = dx_du; - auto dz_dr = -dx_du; + const auto x = mol[i]; + const auto y = mol[i + 1]; + const auto z = mol[i + 2]; + const auto dx_dq = coef * 2.0 * (q * x + u * y - s * z); + const auto dx_dr = coef * 2.0 * (r * x + s * y + u * z); + const auto dy_dr = coef * 2.0 * (s * x - r * y + q * z); + const auto dx_du = coef * 2.0 * (-u * x + q * y + r * z); + const auto dz_ds = dx_dq; + const auto dy_du = -dx_dq; + const auto dy_ds = dx_dr; + const auto dz_du = dx_dr; + const auto dx_ds = -dy_dr; + const auto dz_dq = dy_dr; + const auto dy_dq = dx_du; + const auto dz_dr = -dx_du; gradConverters[k] = dx_dq; gradConverters[k + 1] = dy_dq; gradConverters[k + 2] = dz_dq; @@ -177,17 +178,17 @@ void cartToQuatGrads(const double *quat, const double *mol, int numBPts, } // namespace // atoms/shape features -double calcVolAndGrads(const double *ref, int numRefPts, +double calcVolAndGrads(const double *ref, const int numRefPts, const boost::dynamic_bitset<> *refCarbonRadii, - const double *fit, int numFitPts, + const double *fit, const int numFitPts, const boost::dynamic_bitset<> *fitCarbonRadii, - std::vector &gradConverters, bool useCutoff, - double distCutoff2, const double *quat, - double *gradients) { + std::vector &gradConverters, + const bool useCutoff, const double distCutoff2, + const double *quat, double *gradients) { if (gradients) { cartToQuatGrads(quat, fit, numFitPts, gradConverters, 0); } - static const double CARBON_A = KAPPA / (1.7 * 1.7); + static constexpr double CARBON_A = KAPPA / (1.7 * 1.7); static const double CARBON_BIT = 8.0 * pow(std::numbers::pi / (2 * CARBON_A), 1.5); double vol = 0.0; @@ -195,31 +196,30 @@ double calcVolAndGrads(const double *ref, int numRefPts, // If either of the carbon radii flags aren't supplied, treat them // both as being all carbon. There isn't enough information to do // otherwise. - bool allCarbon = !refCarbonRadii || !fitCarbonRadii; + const bool allCarbon = !refCarbonRadii || !fitCarbonRadii; for (int i = 0, i_idx = 0; i < numRefPts * 4; i += 4, i_idx++) { const auto ai = ref[i + 3]; for (int j = 0, j_idx = 0, k = 0; j < numFitPts * 4; j += 4, j_idx++, k += 12) { - auto dx = ref[i] - fit[j]; - auto dy = ref[i + 1] - fit[j + 1]; - auto dz = ref[i + 2] - fit[j + 2]; - auto d2 = dx * dx + dy * dy + dz * dz; + const auto dx = ref[i] - fit[j]; + const auto dy = ref[i + 1] - fit[j + 1]; + const auto dz = ref[i + 2] - fit[j + 2]; + const auto d2 = dx * dx + dy * dy + dz * dz; if (useCutoff && d2 > distCutoff2) { continue; } const auto aj = fit[j + 3]; - auto mult = -(ai * aj) / (ai + aj); - auto kij = exp(mult * d2); - if (allCarbon || (!allCarbon && (*refCarbonRadii)[i_idx] && - (*fitCarbonRadii)[j_idx])) { + const auto mult = -(ai * aj) / (ai + aj); + const auto kij = exp(mult * d2); + if (allCarbon || ((*refCarbonRadii)[i_idx] && (*fitCarbonRadii)[j_idx])) { vij = kij * CARBON_BIT; } else { - auto pi_ai_aj = std::numbers::pi / (ai + aj); + const auto pi_ai_aj = std::numbers::pi / (ai + aj); vij = 8 * kij * pi_ai_aj * std::sqrt(pi_ai_aj); } vol += vij; if (gradients) { - auto r = 2.0 * vij * mult; + const auto r = 2.0 * vij * mult; // Use the gradient converters to calculate the gradients in quaternion // space. // The zeroth gradient is never used, so don't waste time calculating @@ -247,9 +247,11 @@ double calcVolAndGrads(const double *ref, int numRefPts, } // color features -double calcVolAndGrads(const double *ref, int numRefPts, const int *refTypes, - const double *fit, int numFitPts, const int *fitTypes, - int numFitShape, std::vector &gradConverters, +double calcVolAndGrads(const double *ref, const int numRefPts, + const int *refTypes, const double *fit, + const int numFitPts, const int *fitTypes, + const int numFitShape, + std::vector &gradConverters, const bool useCutoff, const double distCutoff2, const double *quat, double *gradients) { double vol = 0.0; @@ -266,22 +268,22 @@ double calcVolAndGrads(const double *ref, int numRefPts, const int *refTypes, if (aType != bType) { continue; } - auto dx = ref[i] - fit[j]; - auto dy = ref[i + 1] - fit[j + 1]; - auto dz = ref[i + 2] - fit[j + 2]; - auto d2 = dx * dx + dy * dy + dz * dz; + const auto dx = ref[i] - fit[j]; + const auto dy = ref[i + 1] - fit[j + 1]; + const auto dz = ref[i + 2] - fit[j + 2]; + const auto d2 = dx * dx + dy * dy + dz * dz; if (useCutoff && d2 > distCutoff2) { continue; } const auto aj = fit[j + 3]; - auto mult = -(ai * aj) / (ai + aj); - auto kij = exp(mult * d2); + const auto mult = -(ai * aj) / (ai + aj); + const auto kij = exp(mult * d2); - auto pi_ai_aj = std::numbers::pi / (ai + aj); - auto vij = 8 * kij * pi_ai_aj * std::sqrt(pi_ai_aj); + const auto pi_ai_aj = std::numbers::pi / (ai + aj); + const auto vij = 8 * kij * pi_ai_aj * std::sqrt(pi_ai_aj); vol += vij; if (gradients) { - auto r = 2.0 * vij * mult; + const auto r = 2.0 * vij * mult; // Use the converters to calculate the gradients in quaternion space. // The zeroth gradient is never used, so don't waste time calculating // it but leave the code here for completeness and possible future use. @@ -310,7 +312,7 @@ double calcVolAndGrads(const double *ref, int numRefPts, const int *refTypes, void SingleConformerAlignment::applyQuatTrans( const std::array &quatTrans) { // Leave fit at the origin, and move ref to meet it. - RDGeom::Point3D translateA{-quatTrans[4], -quatTrans[5], -quatTrans[6]}; + RDGeom::Point3D const translateA{-quatTrans[4], -quatTrans[5], -quatTrans[6]}; translateShape(d_ref.data(), d_refTemp.data(), d_nRefShape + d_nRefColor, translateA); // Rotate fit by quaternion @@ -343,13 +345,13 @@ void SingleConformerAlignment::calcVolumeAndGradients( d_useCutoff, d_distCutoff2, quatTrans.data(), colorGrads.data()); // The color gradients are normally dwarfed by the shape gradients, so // normalize them and then mix by the same rule as the final score. - auto shapeSum = sqrt(std::accumulate( + const auto shapeSum = sqrt(std::accumulate( gradients.begin() + 1, gradients.end(), 0.0, [](const auto init, const auto g) -> double { return init + g * g; })); - auto colorSum = sqrt(std::accumulate( + const auto colorSum = sqrt(std::accumulate( colorGrads.begin() + 1, colorGrads.end(), 0.0, [](const auto init, const auto g) -> double { return init + g * g; })); - auto ratio = shapeSum / colorSum; + const auto ratio = shapeSum / colorSum; std::transform( gradients.begin() + 1, gradients.end(), colorGrads.begin(), gradients.begin() + 1, [&](const auto g1, const auto g2) -> double { @@ -361,9 +363,9 @@ void SingleConformerAlignment::calcVolumeAndGradients( } bool SingleConformerAlignment::doOverlay(std::array &scores, - unsigned int cycle) { - unsigned int maxIters = cycle == 0 ? 10 : d_maxIts - 10; - auto res = optimise(maxIters); + const unsigned int cycle) { + const unsigned int maxIters = cycle == 0 ? 10 : d_maxIts - 10; + const auto res = optimise(maxIters); // Get the final coords for fit into d_fitTemp, and compute the scores RDGeom::Transform3D xform; @@ -373,7 +375,7 @@ bool SingleConformerAlignment::doOverlay(std::array &scores, applyTransformToShape(d_fit.data(), d_fitTemp.data(), d_nFitShape + d_nFitColor, xform); - auto tscores = calcScores(d_ref.data(), d_fitTemp.data(), true); + const auto tscores = calcScores(d_ref.data(), d_fitTemp.data(), true); scores[0] = tscores[0]; scores[1] = tscores[1]; scores[2] = tscores[2]; @@ -398,14 +400,14 @@ bool SingleConformerAlignment::doOverlay(std::array &scores, } namespace { -double oneStep(double grad, double stepSize, double quatTrans, double oldGrad, - double oldQuatTrans) { +double oneStep(const double grad, const double stepSize, const double quatTrans, + const double oldGrad, const double oldQuatTrans) { double step = 0.0; if (std::signbit(grad) != std::signbit(oldGrad)) { - step = (((quatTrans * fabs(oldGrad)) + (oldQuatTrans * fabs(grad))) / - (fabs(oldGrad) + fabs(grad) + fabs(grad))) - + step = (quatTrans * fabs(oldGrad) + oldQuatTrans * fabs(grad)) / + (fabs(oldGrad) + fabs(grad) + fabs(grad)) - quatTrans; - double newStep = stepSize * grad; + const double newStep = stepSize * grad; if (fabs(step) > fabs(newStep)) { // This is definitely what the PubChem code says! I read it as keeping // the sign of step, but the value of newStep. @@ -417,10 +419,11 @@ double oneStep(double grad, double stepSize, double quatTrans, double oldGrad, return step; } -void calcStep(std::array &grad, double qStepSize, double tStepSize, - std::array &oldGrad, std::array &quatTrans, - std::array &oldQuatTrans, unsigned int iter, - std::array &step) { +void calcStep(const std::array &grad, const double qStepSize, + const double tStepSize, const std::array &oldGrad, + const std::array &quatTrans, + const std::array &oldQuatTrans, + const unsigned int iter, std::array &step) { step[0] = 0.0; if (iter == 0) { // 1st iteration, use default step sizes @@ -446,25 +449,25 @@ void calcStep(std::array &grad, double qStepSize, double tStepSize, } } -double constrainStep(double maxStep, double *step, bool checkSize) { - double mStep = std::max({fabs(step[0]), fabs(step[1]), fabs(step[2])}); +double constrainStep(const double maxStep, double *step, const bool checkSize) { + const double mStep = std::max({fabs(step[0]), fabs(step[1]), fabs(step[2])}); if (mStep > maxStep) { - double scaleFactor = maxStep / mStep; - if (fabs(step[0] > maxStep)) { + const double scaleFactor = maxStep / mStep; + if (fabs(step[0]) > maxStep) { step[0] *= scaleFactor; } - if (fabs(step[1] > maxStep)) { + if (fabs(step[1]) > maxStep) { step[1] *= scaleFactor; } - if (fabs(step[2] > maxStep)) { + if (fabs(step[2]) > maxStep) { step[2] *= scaleFactor; } } if (checkSize) { - double quatSquared = + const double quatSquared = step[0] * step[0] + step[1] * step[1] + step[2] * step[2]; if (quatSquared > 1.0) { - double scaleFactor = 1.0 / (2.0 * quatSquared); + const double scaleFactor = 1.0 / (2.0 * quatSquared); step[0] *= scaleFactor; step[1] *= scaleFactor; step[2] *= scaleFactor; @@ -490,13 +493,14 @@ std::array combineQuatTrans(const std::array &q1, return res; } -double oneReduceStep(double grad, double oldGrad, double quatTrans, - double oldQuatTrans, double stepSize, double step) { +double oneReduceStep(const double grad, const double oldGrad, + const double quatTrans, const double oldQuatTrans, + const double stepSize, double step) { if (std::signbit(grad) != std::signbit(oldGrad)) { - step = (((quatTrans * fabs(oldGrad)) + (oldQuatTrans * fabs(grad))) / - (fabs(oldGrad + fabs(grad)))) - + step = (quatTrans * fabs(oldGrad) + oldQuatTrans * fabs(grad)) / + (fabs(oldGrad) + fabs(grad)) - quatTrans; - double newStep = stepSize * grad; + const double newStep = stepSize * grad; if (fabs(step) > fabs(newStep)) { step *= fabs(newStep / step); } @@ -515,11 +519,12 @@ double oneReduceStep(double grad, double oldGrad, double quatTrans, return step; } -void reduceStep(std::array &grad, std::array &oldGrad, - std::array quatTrans, - std::array &oldQuatTrans, unsigned int lineIter, - std::array &step, double &qStepSize, - double &tStepSize) { +void reduceStep(const std::array &grad, + const std::array &oldGrad, + const std::array &quatTrans, + const std::array &oldQuatTrans, + const unsigned int lineIter, std::array &step, + double &qStepSize, double &tStepSize) { if (lineIter == 2) { qStepSize *= 0.1; tStepSize *= 0.1; @@ -579,7 +584,7 @@ bool SingleConformerAlignment::optimise(unsigned int maxIters) { iter, step); // In case we have to backtrack - double oldComboScore = comboScore; + const double oldComboScore = comboScore; oldQuatTrans = d_quatTrans; oldGrad = grad; @@ -588,8 +593,10 @@ bool SingleConformerAlignment::optimise(unsigned int maxIters) { for (unsigned int lineIter = 0; !converged; lineIter++) { // Check that the absolute max step size does not go beyond some // reasonable size - double mqStep = constrainStep(maxQuaternionStep, step.data() + 1, true); - double mtStep = constrainStep(maxTranslationStep, step.data() + 4, false); + const double mqStep = + constrainStep(maxQuaternionStep, step.data() + 1, true); + const double mtStep = + constrainStep(maxTranslationStep, step.data() + 4, false); if (mqStep <= minQuaternionStep && mtStep <= minTranslationStep) { converged = true; comboScore = 0.0; // Make sure we return to the old one. @@ -597,14 +604,14 @@ bool SingleConformerAlignment::optimise(unsigned int maxIters) { } // Calculate the 0th component of the quaternion. Obviously it // relies on the other 3 components being small - double quatSquared = + const double quatSquared = step[1] * step[1] + step[2] * step[2] + step[3] * step[3]; step[0] = sqrt(1.0 - quatSquared); // Update the quaternion with the step, multiplying them. auto newQuatTrans = combineQuatTrans(d_quatTrans, step); calcVolumeAndGradients(newQuatTrans, shapeOvlpVol, colorOvlpVol, grad); - auto scores = calcScores(shapeOvlpVol, colorOvlpVol); - comboScore = scores[0]; + auto newScores = calcScores(shapeOvlpVol, colorOvlpVol); + comboScore = newScores[0]; // if we made a good step, keep the quaternion and we're done if (comboScore > oldComboScore) { d_quatTrans = newQuatTrans; diff --git a/Code/GraphMol/GaussianShape/catch_tests.cpp b/Code/GraphMol/GaussianShape/catch_tests.cpp index 827507213..1dd0af5db 100644 --- a/Code/GraphMol/GaussianShape/catch_tests.cpp +++ b/Code/GraphMol/GaussianShape/catch_tests.cpp @@ -161,7 +161,7 @@ TEST_CASE("basic alignment") { } else { CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.503, 0.005)); CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(0.761, 0.005)); - CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.245, 0.005)); + CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.233, 0.005)); } // Check that a re-score gives the same answer. const auto rescores = GaussianShape::ScoreMolecule( @@ -294,9 +294,9 @@ TEST_CASE("handling molecules with Hs") { RWMol cp(*probe); RDGeom::Transform3D xform; auto scores = GaussianShape::AlignMolecule(*ref, cp); - CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.700, 0.005)); + CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.736, 0.005)); CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(0.834, 0.005)); - CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.566, 0.005)); + CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.635, 0.005)); for (auto i = 0u; i < cp.getNumAtoms(); ++i) { // the failure mode here was that Hs had HUGE coordinates auto pos = cp.getConformer().getAtomPos(i); @@ -440,19 +440,18 @@ TEST_CASE("Iressa onto Tagrisso") { shapeOpts.allCarbonRadii = false; auto scores = GaussianShape::AlignMolecule(*tagrisso, *iressa, shapeOpts, shapeOpts, nullptr, opts); - CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.332, 0.005)); - CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(0.569, 0.005)); - CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.095, 0.005)); - + CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.331, 0.005)); + CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(0.572, 0.005)); + CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.090, 0.005)); auto rescores = GaussianShape::ScoreMolecule(*tagrisso, *iressa, shapeOpts, shapeOpts, opts); CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(rescores[0], 0.005)); CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(rescores[1], 0.005)); CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(rescores[2], 0.005)); auto aligned_iressa = - "COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCOCC1 |(3.34206,-4.82098,0.224565;2.562,-4.29938,-0.849374;1.40029,-3.66768,-0.520786;0.217637,-4.40844,-0.435429;-0.995315,-3.80966,-0.103538;-2.12266,-4.53841,-0.026421;-3.25097,-3.87666,0.301586;-3.3749,-2.56477,0.559961;-2.22774,-1.8651,0.473675;-2.30832,-0.475114,0.738276;-3.33657,0.464391,0.56286;-4.2686,0.303124,-0.466868;-5.2907,1.23626,-0.641364;-5.38531,2.33529,0.212458;-6.37287,3.22379,0.031359;-4.45782,2.50088,1.24127;-4.5695,3.85508,2.29834;-3.43604,1.56746,1.41601;-0.997368,-2.42737,0.145328;0.187601,-1.67548,0.061391;1.37756,-2.30488,-0.27166;2.52893,-1.56544,-0.352965;2.83995,-1.00034,-1.61907;3.59724,0.302963,-1.40382;2.76275,1.31266,-0.622178;1.52096,1.60458,-1.33518;0.667075,2.50553,-0.54864;-0.616491,2.80465,-1.31789;-0.312021,3.38752,-2.58555;0.491386,2.50505,-3.3701;1.80144,2.19743,-2.65039)|"_smiles; + "COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCOCC1 |(3.08453,-4.88523,0.170077;2.29914,-4.34091,-0.888589;1.15638,-3.68668,-0.538328;-0.040292,-4.40272,-0.438557;-1.23492,-3.78034,-0.0843727;-2.37618,-4.48548,0.00635958;-3.48492,-3.80191,0.355683;-3.57688,-2.48924,0.623248;-2.41657,-1.81353,0.522875;-2.46344,-0.42373,0.79642;-3.47407,0.538362,0.642364;-4.42522,0.402938,-0.373507;-5.42983,1.35852,-0.526776;-5.48774,2.45427,0.334513;-6.45888,3.36454,0.173853;-4.54106,2.59412,1.34956;-4.60745,3.94414,2.41576;-3.53677,1.63825,1.50308;-1.2037,-2.3998,0.172165;-0.00444866,-1.6727,0.0738388;1.16652,-2.3252,-0.281314;2.33196,-1.60985,-0.376552;2.6352,-1.04403,-1.64422;3.4233,0.241654,-1.43364;2.62275,1.26418,-0.633414;1.37651,1.58654,-1.32523;0.554317,2.50073,-0.520411;-0.734405,2.83149,-1.26777;-0.437359,3.41523,-2.53679;0.334793,2.52057,-3.33873;1.64904,2.18104,-2.64136)|"_smiles; REQUIRE(aligned_iressa); - checkMolsHaveRoughlySameCoords(*iressa, *aligned_iressa); + CHECK(checkMolsHaveRoughlySameCoords(*iressa, *aligned_iressa)); } TEST_CASE("Optimise in place") { @@ -501,9 +500,9 @@ TEST_CASE("Optimise in place") { ROMol cp(*canon_probe); auto scores = GaussianShape::AlignMolecule(*pdb_trp_3tmn, cp, shapeOpts, shapeOpts, nullptr, opts); - CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.197, 0.001)); - CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(0.361, 0.001)); - CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.033, 0.001)); + CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.201, 0.001)); + CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(0.364, 0.001)); + CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.039, 0.001)); } { // And with reference as a shape the same @@ -561,10 +560,9 @@ TEST_CASE("Fragment Mode") { // Use the smaller molecule as the probe auto scores = GaussianShape::AlignShape(refShape, probeShape, &xform, opts); // These are close to the values above for starting from the xtal structures. - CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.311, 0.005)); - CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(0.408, 0.005)); - CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.215, 0.005)); - MolTransforms::transformConformer(pdb_trp_3tmn->getConformer(), xform); + CHECK_THAT(scores[0], Catch::Matchers::WithinAbs(0.332, 0.005)); + CHECK_THAT(scores[1], Catch::Matchers::WithinAbs(0.413, 0.005)); + CHECK_THAT(scores[2], Catch::Matchers::WithinAbs(0.251, 0.005)); } TEST_CASE("custom feature points") {