mirror of
https://github.com/rdkit/rdkit.git
synced 2026-06-03 21:44:30 +08:00
9f551aedbe
* First import of GaussianShape. * Tidying. * Custom features. * Optimise. * Optimise. * Return 3 scores rather than 2 including combo score. * Rename useFeatures to useColors. * Python wrappers. * Python tests. * Take out big test. * Add new start mode, as PubChem does it. * Doh! * Fix MolTransforms eigenvalue return. * Two cycle optimisation, mostly working. * Take out bestSoFar score from SCA. * Take out DTYPE. * Tidy out redundant variables. * Optimisation in 2 parts. * More fiddling in pursuit of speed. * Update Python wrapper. * Tweak. * Atom subsets and different radii. * Fix test. * Revert pubchem_shape's test.cpp. * Serialize ShapeInput. * Trigger build * Remove pointers to std::arrays in ShapeInput. * ShapeInput virtual d'tor. * Precondition - ShapeInput needs a molecule with at least 1 conformer. * Rename ShapeInput::d_centroid to ShapeInput::d_canonTrans. * Fix normalization bugs. * Select start mode using moments of inertia rather than eigenvalues of canonical transformation. * Include color features in moments of inertia. * Smidge faster. * Tversky similarity. * Tidy tests. * Tests working on Linux. * Revert force of right handed axes in MolTransforms::computePrincipalAxesAndMomentsFromGyrationMatrix replacing with a comment in the code. * Response to review. * Sneaky allCarbon bug. * add multithreaded test * Response to review. * Doh! Don't recalculate normalization after every transformation. * Re-instate d_normalizationOK. * Re-name functions for fetching canonical transformations. * Separate alpha from coords. * MultiConf works with single conf extraction. * Extract all conformations. Max and best similarities. * Renames d_currConformer to d_activeShape. * Update shapeToMol. * Update shapeToMol. * Changes from synthon shape searching. * Fix normalization of multiple confs. * Update Python wrappers. * Fix shape merge. * Improve bestSimilarity. * Fix python wrapper. * Pull in changes from SynthonShapeSearch: make pruneShapes public. function to negate Alpha values. * clang-tidy suggestions. * clang-tidy suggestions. * Bug in quaternion gradients - we now have only 3 coordinates. * Tidy tests. * Mac result slightly different. * Multi conformer molecule alignment. * Optionally return raw overlap volumes in score functions. * Python wrappers for raw overlap volumes. * Update Python wrapper ShapeInputOptions. * Tidy for PR. * Extra include file. * Extra library * Tidy forward declarations. * Don't prune if threshold < 0.0. * Windows exporty thing. * Check SMILES on merge of ShapeInputs. * PRECONDITION of SMILES on merge of ShapeInputs. * Response to review - rename some functions. * change how overlapVols is passed add a test for it * API suggestions * Response to review. * Remove debugging writes. * Fix Python wrappers. --------- Co-authored-by: David Cosgrove <david@cozchemix.co.uk> Co-authored-by: greg landrum <greg.landrum@gmail.com>
543 lines
23 KiB
C++
543 lines
23 KiB
C++
//
|
|
// Copyright (C) 2026 David Cosgrove 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.
|
|
//
|
|
// Original author: David Cosgrove (CozChemIx Limited)
|
|
//
|
|
// This is an implementation of the Gaussian overlap molecular overlay
|
|
// method of Grant, Pickup and Gallardo.
|
|
// J. Comp. Chem., 17, 1653-1666 (1996)
|
|
// https://doi.org/10.1002/(SICI)1096-987X(19961115)17:14%3C1653::AID-JCC7%3E3.0.CO;2-K
|
|
// It uses implementation ideas and some code from the PubChem implementation
|
|
// https://github.com/ncbi/pubchem-align3d/blob/main/shape_neighbor.cpp.
|
|
|
|
#include <cmath>
|
|
#include <numbers>
|
|
|
|
#include <Geometry/Transform3D.h>
|
|
#include <GraphMol/ROMol.h>
|
|
#include <GraphMol/GaussianShape/GaussianShape.h>
|
|
|
|
#include "GraphMol/SmilesParse/SmilesWrite.h"
|
|
|
|
#include <GraphMol/GaussianShape/ShapeInput.h>
|
|
#include <GraphMol/GaussianShape/SingleConformerAlignment.h>
|
|
|
|
namespace RDKit {
|
|
namespace GaussianShape {
|
|
|
|
namespace {
|
|
// Compute final overlay transform, which applies fitShape's
|
|
// initial canonical transformation, followed by the overlay transform and
|
|
// finally the inverse of refShape's initial canonical transformation.
|
|
RDGeom::Transform3D computeFinalTransform(
|
|
const std::array<double, 3> &inRefTrans,
|
|
const std::array<double, 9> &inRefRot,
|
|
const std::array<double, 3> &inFitTrans,
|
|
const std::array<double, 9> &inFitRot, const RDGeom::Transform3D &ovXform) {
|
|
// Move to fitShape's initial centroid and principal axes
|
|
RDGeom::Transform3D transform0;
|
|
transform0.SetTranslation(
|
|
RDGeom::Point3D{inFitTrans[0], inFitTrans[1], inFitTrans[2]});
|
|
|
|
RDGeom::Transform3D transform1;
|
|
transform1.setValUnchecked(0, 0, inFitRot[0]);
|
|
transform1.setValUnchecked(0, 1, inFitRot[1]);
|
|
transform1.setValUnchecked(0, 2, inFitRot[2]);
|
|
transform1.setValUnchecked(1, 0, inFitRot[3]);
|
|
transform1.setValUnchecked(1, 1, inFitRot[4]);
|
|
transform1.setValUnchecked(1, 2, inFitRot[5]);
|
|
transform1.setValUnchecked(2, 0, inFitRot[6]);
|
|
transform1.setValUnchecked(2, 1, inFitRot[7]);
|
|
transform1.setValUnchecked(2, 2, inFitRot[8]);
|
|
|
|
RDGeom::Transform3D toRefRefFrame;
|
|
// Rotate by the inverse of the ref shape's canonical rotation and
|
|
// translate by the negative of its canonical translation.
|
|
toRefRefFrame.setValUnchecked(0, 0, inRefRot[0]);
|
|
toRefRefFrame.setValUnchecked(0, 1, inRefRot[3]);
|
|
toRefRefFrame.setValUnchecked(0, 2, inRefRot[6]);
|
|
toRefRefFrame.setValUnchecked(0, 3, -inRefTrans[0]);
|
|
toRefRefFrame.setValUnchecked(1, 0, inRefRot[1]);
|
|
toRefRefFrame.setValUnchecked(1, 1, inRefRot[4]);
|
|
toRefRefFrame.setValUnchecked(1, 2, inRefRot[7]);
|
|
toRefRefFrame.setValUnchecked(1, 3, -inRefTrans[1]);
|
|
toRefRefFrame.setValUnchecked(2, 0, inRefRot[2]);
|
|
toRefRefFrame.setValUnchecked(2, 1, inRefRot[5]);
|
|
toRefRefFrame.setValUnchecked(2, 2, inRefRot[8]);
|
|
toRefRefFrame.setValUnchecked(2, 3, -inRefTrans[2]);
|
|
|
|
auto finalTransform = toRefRefFrame * ovXform * transform1 * transform0;
|
|
return finalTransform;
|
|
}
|
|
|
|
// Return the original transformation quaternion for the given index.
|
|
// Different optimisation modes have different numbers of starting
|
|
// orientations to try. In order these are no transformation, rotate 180
|
|
// degrees about each axis and rotate +/- 45 degrees about 2 axes at a time.
|
|
std::array<double, 4> getInitialRotationPlain(
|
|
int index, const ShapeInput &refShape, const ShapeInput &fitShape,
|
|
const RDGeom::Point3D &refDisp, const ShapeOverlayOptions &overlayOpts,
|
|
double &score) {
|
|
static const double sinpi_4 = std::sin(std::numbers::pi / 4.0);
|
|
const static std::vector<std::array<double, 4>> quats{
|
|
{1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0},
|
|
{0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 1.0},
|
|
{sinpi_4, -sinpi_4, 0.0, 0.0}, {sinpi_4, sinpi_4, 0.0, 0.0},
|
|
{0.0, 0.0, -sinpi_4, sinpi_4}, {0.0, 0.0, sinpi_4, sinpi_4},
|
|
{sinpi_4, 0.0, 0.0, -sinpi_4}, {0.0, sinpi_4, sinpi_4, 0.0},
|
|
{sinpi_4, 0.0, 0.0, sinpi_4}, {0.0, -sinpi_4, sinpi_4, 0.0},
|
|
{sinpi_4, 0.0, sinpi_4, 0.0}, {0.0, sinpi_4, 0.0, sinpi_4},
|
|
{0.0, -sinpi_4, 0.0, sinpi_4}, {sinpi_4, 0.0, -sinpi_4, 0.0}};
|
|
const bool useColor = overlayOpts.optimMode != OptimMode::SHAPE_ONLY;
|
|
const std::array<double, 7> quatTrans{
|
|
quats[index][0], quats[index][1], quats[index][2], quats[index][3],
|
|
refDisp[0], refDisp[1], refDisp[2]};
|
|
SingleConformerAlignment sca(
|
|
refShape.getCoords(), refShape.getAlphas(),
|
|
refShape.getFeatureTypes().data(), refShape.getCarbonRadii(),
|
|
refShape.getNumAtoms(), refShape.getNumFeatures(),
|
|
refShape.getShapeVolume(), refShape.getColorVolume(),
|
|
fitShape.getCoords(), fitShape.getAlphas(),
|
|
fitShape.getFeatureTypes().data(), fitShape.getCarbonRadii(),
|
|
fitShape.getNumAtoms(), fitShape.getNumFeatures(),
|
|
fitShape.getShapeVolume(), fitShape.getColorVolume(), quatTrans,
|
|
overlayOpts.optimMode, overlayOpts.simAlpha, overlayOpts.simBeta,
|
|
overlayOpts.optParam, overlayOpts.useDistCutoff, overlayOpts.distCutoff,
|
|
overlayOpts.shapeConvergenceCriterion, overlayOpts.nSteps);
|
|
const auto scores = sca.calcScores(useColor);
|
|
score = scores[0];
|
|
return quats[index];
|
|
}
|
|
|
|
// Return the initial transformation matrix in the manner of the PubChem
|
|
// overlay code. Rotate 180 degrees about each axis, and then
|
|
// add +/ ~25 degrees from that. It is not revealed where that
|
|
// angle comes from.
|
|
std::array<double, 4> getInitialRotationWiggle(
|
|
int index, const ShapeInput &refShape, const ShapeInput &fitShape,
|
|
const RDGeom::Point3D &refDisp, const ShapeOverlayOptions &overlayOpts,
|
|
double &score) {
|
|
const static double qrot1 = 0.977659114061,
|
|
qrot = 0.210196709523; // 0.215 (un-normalized)
|
|
const static std::vector<std::array<double, 4>> quats{
|
|
{1.0, 0.0, 0.0, 0.0}, // 0 X, Y, Z
|
|
{qrot1, qrot, 0.0, 0.0}, {qrot1, -qrot, 0.0, 0.0},
|
|
{qrot1, 0.0, qrot, 0.0}, {qrot1, 0.0, -qrot, 0.0},
|
|
{qrot1, 0.0, 0.0, qrot}, {qrot1, 0.0, 0.0, -qrot},
|
|
{0.0, 1.0, 0.0, 0.0}, // 1 X, -Y, -Z
|
|
{qrot, qrot1, 0.0, 0.0}, {qrot, -qrot1, 0.0, 0.0},
|
|
{0.0, qrot1, qrot, 0.0}, {0.0, qrot1, -qrot, 0.0},
|
|
{0.0, qrot1, 0.0, qrot}, {0.0, qrot1, 0.0, -qrot},
|
|
{0.0, 0.0, 0.0, 1.0}, // 2 -X, -Y, Z
|
|
{qrot, 0.0, 0.0, qrot1}, {qrot, 0.0, 0.0, -qrot1},
|
|
{0.0, qrot, 0.0, qrot1}, {0.0, -qrot, 0.0, qrot1},
|
|
{0.0, 0.0, qrot, qrot1}, {0.0, 0.0, -qrot, qrot1},
|
|
{0.0, 0.0, 1.0, 0.0}, // 3 -X, Y, -Z
|
|
{qrot, 0.0, qrot1, 0.0}, {qrot, 0.0, -qrot1, 0.0},
|
|
{0.0, qrot, qrot1, 0.0}, {0.0, -qrot, qrot1, 0.0},
|
|
{0.0, 0.0, qrot1, qrot}, {0.0, 0.0, qrot1, -qrot}};
|
|
unsigned int start_quat = index * 7;
|
|
unsigned int bestQuat = 0;
|
|
double bestScore = 0.0;
|
|
bool useColor = overlayOpts.optimMode != OptimMode::SHAPE_ONLY;
|
|
std::array<double, 7> tmpQuatTrans{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
|
|
SingleConformerAlignment sca(
|
|
refShape.getCoords(), refShape.getAlphas(),
|
|
refShape.getFeatureTypes().data(), refShape.getCarbonRadii(),
|
|
refShape.getNumAtoms(), refShape.getNumFeatures(),
|
|
refShape.getShapeVolume(), refShape.getColorVolume(),
|
|
fitShape.getCoords(), fitShape.getAlphas(),
|
|
fitShape.getFeatureTypes().data(), fitShape.getCarbonRadii(),
|
|
fitShape.getNumAtoms(), fitShape.getNumFeatures(),
|
|
fitShape.getShapeVolume(), fitShape.getColorVolume(), tmpQuatTrans,
|
|
overlayOpts.optimMode, overlayOpts.simAlpha, overlayOpts.simBeta,
|
|
overlayOpts.optParam, overlayOpts.useDistCutoff, overlayOpts.distCutoff,
|
|
overlayOpts.shapeConvergenceCriterion, overlayOpts.nSteps);
|
|
|
|
for (unsigned int i = start_quat; i < start_quat + 7; ++i) {
|
|
std::array<double, 7> quatTrans{quats[i][0], quats[i][1], quats[i][2],
|
|
quats[i][3], refDisp[0], refDisp[1],
|
|
refDisp[2]};
|
|
sca.setQuatTrans(quatTrans);
|
|
auto scores = sca.calcScores(useColor);
|
|
if (scores[0] > bestScore) {
|
|
bestScore = scores[0];
|
|
bestQuat = i;
|
|
}
|
|
}
|
|
score = bestScore;
|
|
return quats[bestQuat];
|
|
}
|
|
|
|
// Return the translation that puts the extreme of refShape at the
|
|
// extreme of the fitShape along the appropriate axis.
|
|
RDGeom::Point3D getInitialTranslation(int index, ShapeInput &refShape,
|
|
ShapeInput fitShape) {
|
|
auto getDisp = [](ShapeInput &shape, size_t i) -> RDGeom::Point3D {
|
|
const double *coord =
|
|
shape.getCoords().data() + shape.calcExtremes()[i] * 3;
|
|
return RDGeom::Point3D(coord[0], coord[1], coord[2]);
|
|
};
|
|
RDGeom::Point3D disp;
|
|
RDGeom::Point3D refDisp, fitDisp;
|
|
switch (index) {
|
|
case 1:
|
|
refDisp = getDisp(refShape, 0);
|
|
fitDisp = getDisp(fitShape, 0);
|
|
disp = fitDisp - refDisp;
|
|
break;
|
|
case 2:
|
|
refDisp = getDisp(refShape, 1);
|
|
fitDisp = getDisp(fitShape, 1);
|
|
disp = fitDisp - refDisp;
|
|
break;
|
|
case 3:
|
|
refDisp = getDisp(refShape, 2);
|
|
fitDisp = getDisp(fitShape, 2);
|
|
disp = fitDisp - refDisp;
|
|
break;
|
|
case 4:
|
|
refDisp = getDisp(refShape, 3);
|
|
fitDisp = getDisp(fitShape, 3);
|
|
disp = fitDisp - refDisp;
|
|
break;
|
|
case 5:
|
|
refDisp = getDisp(refShape, 4);
|
|
fitDisp = getDisp(fitShape, 4);
|
|
disp = fitDisp - refDisp;
|
|
break;
|
|
case 6:
|
|
refDisp = getDisp(refShape, 5);
|
|
fitDisp = getDisp(fitShape, 5);
|
|
disp = fitDisp - refDisp;
|
|
break;
|
|
default:
|
|
break;
|
|
}
|
|
return disp;
|
|
}
|
|
|
|
// This is how the PubChem code decides between ROTATE_180_WIGGLE and
|
|
// ROTATE_45. I have no clue.
|
|
unsigned int calculateQrat(const std::array<double, 3> &eigenValues) {
|
|
double double_ev_oe[3]{eigenValues[1] + eigenValues[2] - eigenValues[0],
|
|
eigenValues[0] + eigenValues[2] - eigenValues[1],
|
|
eigenValues[0] + eigenValues[1] - eigenValues[2]};
|
|
std::sort(double_ev_oe, double_ev_oe + 3, std::greater<double>());
|
|
|
|
constexpr static double qrat_threshold = 0.7225; // 0.85*0.85;
|
|
unsigned int qrat = 1000;
|
|
|
|
if (double_ev_oe[1] > 0) {
|
|
unsigned int u_rqyx, u_rqzy;
|
|
if (qrat_threshold < double_ev_oe[1] / double_ev_oe[0]) {
|
|
u_rqyx = 1;
|
|
} else {
|
|
u_rqyx = 0;
|
|
}
|
|
if (qrat_threshold < double_ev_oe[2] / double_ev_oe[1]) {
|
|
u_rqzy = 1;
|
|
} else {
|
|
u_rqzy = 0;
|
|
}
|
|
|
|
qrat = u_rqyx + u_rqzy;
|
|
}
|
|
return qrat;
|
|
}
|
|
|
|
StartMode decideStartModeFromEigenValues(const ShapeInput &refShape,
|
|
const ShapeInput &fitShape) {
|
|
// The PubChem code uses the moments of inertia for this, rather than the
|
|
// canonical transformation.
|
|
const auto rqratwf = calculateQrat(refShape.calcMomentsOfInertia(true));
|
|
const auto fqratwf = calculateQrat(fitShape.calcMomentsOfInertia(true));
|
|
StartMode startModeWF{StartMode::ROTATE_180_WIGGLE};
|
|
if (rqratwf > 0 || fqratwf > 0) {
|
|
startModeWF = StartMode::ROTATE_45;
|
|
}
|
|
return startModeWF;
|
|
}
|
|
|
|
std::array<double, 3> alignShape(ShapeInput &refShape, ShapeInput &fitShape,
|
|
RDGeom::Transform3D &bestXform,
|
|
const ShapeOverlayOptions &overlayOpts) {
|
|
unsigned int finalRotIndex = 1;
|
|
auto startMode = overlayOpts.startMode;
|
|
if (startMode == StartMode::A_LA_PUBCHEM) {
|
|
startMode = decideStartModeFromEigenValues(refShape, fitShape);
|
|
}
|
|
|
|
switch (startMode) {
|
|
case StartMode::ROTATE_0:
|
|
case StartMode::ROTATE_0_FRAGMENT:
|
|
break;
|
|
case StartMode::ROTATE_180:
|
|
case StartMode::ROTATE_180_FRAGMENT:
|
|
case StartMode::ROTATE_180_WIGGLE:
|
|
finalRotIndex = 4;
|
|
break;
|
|
case StartMode::ROTATE_45:
|
|
case StartMode::ROTATE_45_FRAGMENT:
|
|
finalRotIndex = 16;
|
|
break;
|
|
default:
|
|
break;
|
|
}
|
|
unsigned int finalTransIndex = 1;
|
|
if (startMode == StartMode::ROTATE_0_FRAGMENT ||
|
|
startMode == StartMode::ROTATE_45_FRAGMENT ||
|
|
startMode == StartMode::ROTATE_180_FRAGMENT) {
|
|
finalTransIndex = 7;
|
|
}
|
|
|
|
std::array<double, 3> bestScore;
|
|
double bestTotal = -1.0;
|
|
|
|
// Get together the start transformations.
|
|
std::vector<std::unique_ptr<SingleConformerAlignment>> aligners;
|
|
std::vector<std::pair<double, unsigned int>> bestScoreForStart;
|
|
bestScoreForStart.reserve(finalTransIndex * finalRotIndex);
|
|
unsigned int k = 0;
|
|
for (unsigned int j = 0; j < finalTransIndex; j++) {
|
|
auto refDisp = getInitialTranslation(j, refShape, fitShape);
|
|
std::array<double, 4> quat;
|
|
for (unsigned int i = 0; i < finalRotIndex; i++, k++) {
|
|
double score = 0.0;
|
|
if (startMode == StartMode::ROTATE_180_WIGGLE) {
|
|
quat = getInitialRotationWiggle(i, refShape, fitShape, refDisp,
|
|
overlayOpts, score);
|
|
} else {
|
|
quat = getInitialRotationPlain(i, refShape, fitShape, refDisp,
|
|
overlayOpts, score);
|
|
}
|
|
std::array<double, 7> initQuat{quat[0], quat[1], quat[2], quat[3],
|
|
refDisp.x, refDisp.y, refDisp.z};
|
|
aligners.emplace_back(std::make_unique<SingleConformerAlignment>(
|
|
refShape.getCoords(), refShape.getAlphas(),
|
|
refShape.getFeatureTypes().data(), refShape.getCarbonRadii(),
|
|
refShape.getNumAtoms(), refShape.getNumFeatures(),
|
|
refShape.getShapeVolume(), refShape.getColorVolume(),
|
|
fitShape.getCoords(), fitShape.getAlphas(),
|
|
fitShape.getFeatureTypes().data(), fitShape.getCarbonRadii(),
|
|
fitShape.getNumAtoms(), fitShape.getNumFeatures(),
|
|
fitShape.getShapeVolume(), fitShape.getColorVolume(), initQuat,
|
|
overlayOpts.optimMode, overlayOpts.simAlpha, overlayOpts.simBeta,
|
|
overlayOpts.optParam, overlayOpts.useDistCutoff,
|
|
overlayOpts.distCutoff, overlayOpts.shapeConvergenceCriterion,
|
|
overlayOpts.nSteps));
|
|
bestScoreForStart.push_back({score, k});
|
|
}
|
|
}
|
|
|
|
// Do it in 2 cycles, a quick optimisation first, followed by an additional
|
|
// longer one for those that look like they're going to win.
|
|
for (unsigned int cycle = 0; cycle < 2; cycle++) {
|
|
std::ranges::sort(bestScoreForStart,
|
|
[](const auto &p1, const auto &p2) -> bool {
|
|
return p1.first > p2.first;
|
|
});
|
|
std::vector<std::pair<double, unsigned int>> nextBestScoreForStart;
|
|
nextBestScoreForStart.reserve(finalTransIndex * finalRotIndex);
|
|
for (const auto &[bssf, m] : bestScoreForStart) {
|
|
if (cycle == 1) {
|
|
if (bssf < 0.7 * bestScore[0]) {
|
|
continue;
|
|
}
|
|
}
|
|
std::array<double, 20> outScores;
|
|
aligners[m]->doOverlay(outScores, cycle);
|
|
nextBestScoreForStart.emplace_back(outScores[0], m);
|
|
if (outScores[0] > bestTotal) {
|
|
bestTotal = outScores[0];
|
|
bestScore =
|
|
std::array<double, 3>{outScores[0], outScores[1], outScores[2]};
|
|
aligners[m]->getFinalQuatTrans(bestXform);
|
|
}
|
|
}
|
|
bestScoreForStart = nextBestScoreForStart;
|
|
}
|
|
return bestScore;
|
|
}
|
|
} // namespace
|
|
|
|
std::array<double, 3> AlignShape(const ShapeInput &refShape,
|
|
ShapeInput &fitShape,
|
|
RDGeom::Transform3D *xform,
|
|
const ShapeOverlayOptions &overlayOpts) {
|
|
// The shapes aren't necessarily normalized (it's not done on creation, for
|
|
// example) but they might need to be.
|
|
auto workingRefShape = std::make_unique<ShapeInput>(refShape);
|
|
auto workingFitShape = std::make_unique<ShapeInput>(fitShape);
|
|
const auto inRefTrans = workingRefShape->calcCanonicalTranslation();
|
|
const auto inRefRot = workingRefShape->calcCanonicalRotation();
|
|
const auto inFitTrans = workingFitShape->calcCanonicalTranslation();
|
|
const auto inFitRot = workingFitShape->calcCanonicalRotation();
|
|
// If we're not normalizing, translate both shapes so that the fit
|
|
// is at the origin, so the rotations work.
|
|
RDGeom::Transform3D moveToOrigin;
|
|
RDGeom::Transform3D moveFromOrigin;
|
|
if (overlayOpts.normalize) {
|
|
if (!workingRefShape->getIsNormalized()) {
|
|
workingRefShape->normalizeCoords();
|
|
}
|
|
if (!workingFitShape->getIsNormalized()) {
|
|
workingFitShape->normalizeCoords();
|
|
}
|
|
} else {
|
|
const auto &canonTrans = workingFitShape->calcCanonicalTranslation();
|
|
moveToOrigin.SetTranslation(
|
|
RDGeom::Point3D{canonTrans[0], canonTrans[1], canonTrans[2]});
|
|
moveFromOrigin.SetTranslation(
|
|
RDGeom::Point3D{-canonTrans[0], -canonTrans[1], -canonTrans[2]});
|
|
workingFitShape->transformCoords(moveToOrigin);
|
|
workingRefShape->transformCoords(moveToOrigin);
|
|
}
|
|
|
|
RDGeom::Transform3D bestXform;
|
|
const auto scores =
|
|
alignShape(*workingRefShape, *workingFitShape, bestXform, overlayOpts);
|
|
if (!overlayOpts.normalize) {
|
|
// Shove it back again.
|
|
bestXform = moveFromOrigin * bestXform * moveToOrigin;
|
|
} else {
|
|
bestXform = computeFinalTransform(inRefTrans, inRefRot, inFitTrans,
|
|
inFitRot, bestXform);
|
|
}
|
|
fitShape.transformCoords(bestXform);
|
|
if (xform) {
|
|
*xform = bestXform;
|
|
}
|
|
|
|
return scores;
|
|
}
|
|
|
|
std::array<double, 3> AlignMolecule(const ShapeInput &refShape, ROMol &fit,
|
|
const ShapeInputOptions &fitOpts,
|
|
RDGeom::Transform3D *xform,
|
|
const ShapeOverlayOptions &overlayOpts,
|
|
int fitConfId) {
|
|
auto fitShape = ShapeInput(fit, fitConfId, fitOpts, overlayOpts);
|
|
RDGeom::Transform3D tmpXform;
|
|
const auto scores = AlignShape(refShape, fitShape, &tmpXform, overlayOpts);
|
|
MolTransforms::transformConformer(fit.getConformer(fitConfId), tmpXform);
|
|
if (xform) {
|
|
*xform = tmpXform;
|
|
}
|
|
return scores;
|
|
}
|
|
|
|
std::array<double, 3> AlignMolecule(const ROMol &ref, ROMol &fit,
|
|
const ShapeInputOptions &refOpts,
|
|
const ShapeInputOptions &fitOpts,
|
|
RDGeom::Transform3D *xform,
|
|
const ShapeOverlayOptions &overlayOpts,
|
|
int refConfId, int fitConfId) {
|
|
const auto refShape = ShapeInput(ref, refConfId, refOpts, overlayOpts);
|
|
const auto scores =
|
|
AlignMolecule(refShape, fit, fitOpts, xform, overlayOpts, fitConfId);
|
|
return scores;
|
|
}
|
|
|
|
void ScoreMoleculeAllConformers(const ROMol &ref, const ROMol &fit,
|
|
int &refConfId, int &fitConfId,
|
|
std::vector<std::vector<double>> &combScores,
|
|
const ShapeInputOptions &refOpts,
|
|
const ShapeInputOptions &fitOpts,
|
|
const ShapeOverlayOptions &overlayOpts,
|
|
RDGeom::Transform3D *xform) {
|
|
// Pruning the shapes wastes time and obviously removes the correspondence
|
|
// between conformers and shapes.
|
|
auto refOptsCp = refOpts;
|
|
refOptsCp.shapePruneThreshold = -1;
|
|
refOptsCp.sortShapes = false;
|
|
auto fitOptsCp = fitOpts;
|
|
fitOptsCp.shapePruneThreshold = -1;
|
|
fitOptsCp.sortShapes = false;
|
|
auto refShape = ShapeInput(ref, -1, refOptsCp, overlayOpts);
|
|
auto fitShape = ShapeInput(fit, -1, fitOptsCp, overlayOpts);
|
|
combScores = std::vector<std::vector<double>>(
|
|
refShape.getNumShapes(), std::vector<double>(fitShape.getNumShapes()));
|
|
double bestScore = -1.0;
|
|
for (unsigned int i = 0; i < refShape.getNumShapes(); i++) {
|
|
refShape.setActiveShape(i);
|
|
for (unsigned int j = 0; j < fitShape.getNumShapes(); j++) {
|
|
fitShape.setActiveShape(j);
|
|
RDGeom::Transform3D thisXform;
|
|
auto scores = AlignShape(refShape, fitShape, &thisXform, overlayOpts);
|
|
combScores[i][j] = scores[0];
|
|
if (scores[0] > bestScore) {
|
|
bestScore = scores[0];
|
|
refConfId = i;
|
|
fitConfId = j;
|
|
if (xform) {
|
|
*xform = thisXform;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
std::array<double, 3> ScoreShape(const ShapeInput &refShape,
|
|
const ShapeInput &fitShape,
|
|
const ShapeOverlayOptions &overlayOpts,
|
|
std::pair<double, double> *overlapVols) {
|
|
auto refWorking = refShape.getCoords();
|
|
auto fitWorking = fitShape.getCoords();
|
|
std::array<double, 7> quatTrans{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
|
|
SingleConformerAlignment sca(
|
|
refShape.getCoords(), refShape.getAlphas(),
|
|
refShape.getFeatureTypes().data(), refShape.getCarbonRadii(),
|
|
refShape.getNumAtoms(), refShape.getNumFeatures(),
|
|
refShape.getShapeVolume(), refShape.getColorVolume(),
|
|
fitShape.getCoords(), fitShape.getAlphas(),
|
|
fitShape.getFeatureTypes().data(), fitShape.getCarbonRadii(),
|
|
fitShape.getNumAtoms(), fitShape.getNumFeatures(),
|
|
fitShape.getShapeVolume(), fitShape.getColorVolume(), quatTrans,
|
|
overlayOpts.optimMode, overlayOpts.simAlpha, overlayOpts.simBeta,
|
|
overlayOpts.optParam, overlayOpts.useDistCutoff, overlayOpts.distCutoff,
|
|
overlayOpts.shapeConvergenceCriterion, overlayOpts.nSteps);
|
|
const bool includeColor = overlayOpts.optimMode != OptimMode::SHAPE_ONLY;
|
|
const auto scores = sca.calcScores(refShape.getCoords().data(),
|
|
fitShape.getCoords().data(), includeColor);
|
|
if (overlapVols) {
|
|
(*overlapVols) = std::make_pair(scores[3], scores[4]);
|
|
}
|
|
return std::array{scores[0], scores[1], scores[2]};
|
|
}
|
|
|
|
std::array<double, 3> ScoreMolecule(const ShapeInput &refShape,
|
|
const ROMol &fit,
|
|
const ShapeInputOptions &fitOpts,
|
|
const ShapeOverlayOptions &overlayOpts,
|
|
int fitConfId,
|
|
std::pair<double, double> *overlapVols) {
|
|
const auto fitShape = ShapeInput(fit, fitConfId, fitOpts, overlayOpts);
|
|
return ScoreShape(refShape, fitShape, overlayOpts, overlapVols);
|
|
}
|
|
|
|
std::array<double, 3> ScoreMolecule(const ROMol &ref, const ROMol &fit,
|
|
const ShapeInputOptions &refOpts,
|
|
const ShapeInputOptions &fitOpts,
|
|
const ShapeOverlayOptions &overlayOpts,
|
|
int refConfId, int fitConfId,
|
|
std::pair<double, double> *overlapVols) {
|
|
ShapeOverlayOptions tmpOpts = overlayOpts;
|
|
tmpOpts.normalize = false;
|
|
tmpOpts.startMode = StartMode::ROTATE_0;
|
|
ShapeInputOptions tmpRefOpts = refOpts;
|
|
auto refShape = ShapeInput(ref, refConfId, refOpts, tmpOpts);
|
|
|
|
ShapeInputOptions tmpFitOpts = fitOpts;
|
|
const auto fitShape = ShapeInput(fit, fitConfId, fitOpts, tmpOpts);
|
|
|
|
return ScoreShape(refShape, fitShape, tmpOpts, overlapVols);
|
|
}
|
|
} // namespace GaussianShape
|
|
} // namespace RDKit
|