mirror of
https://github.com/rdkit/rdkit.git
synced 2026-06-04 21:54:27 +08:00
d54d0c5a10
* Different test for possible coords inversion. * Typo. * Comment in test. --------- Co-authored-by: David Cosgrove <david@cozchemix.co.uk>
700 lines
25 KiB
C++
700 lines
25 KiB
C++
//
|
|
// Copyright (C) 2003-2024 Greg Landrum 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.
|
|
//
|
|
#include "MolTransforms.h"
|
|
#include <GraphMol/RDKitBase.h>
|
|
#include <GraphMol/QueryOps.h>
|
|
#include <Numerics/EigenSolvers/PowerEigenSolver.h>
|
|
#include <Numerics/SymmMatrix.h>
|
|
#include <Numerics/Matrix.h>
|
|
#include <Geometry/Transform3D.h>
|
|
#include <stack>
|
|
#include <boost/dynamic_bitset.hpp>
|
|
#include <RDGeneral/Exceptions.h>
|
|
|
|
#ifndef RDK_HAS_EIGEN3
|
|
constexpr double EIGEN_TOLERANCE = 5.0e-2;
|
|
#endif
|
|
namespace MolTransforms {
|
|
|
|
using namespace RDKit;
|
|
void transformAtom(Atom *atom, RDGeom::Transform3D &tform) {
|
|
PRECONDITION(atom, "no atom");
|
|
ROMol &mol = atom->getOwningMol();
|
|
for (ROMol::ConstConformerIterator ci = mol.beginConformers();
|
|
ci != mol.endConformers(); ci++) {
|
|
RDGeom::Point3D &pos = (*ci)->getAtomPos(atom->getIdx());
|
|
tform.TransformPoint(pos);
|
|
}
|
|
// atom->setPos(pos);
|
|
}
|
|
void transformMolsAtoms(ROMol *mol, RDGeom::Transform3D &tform) {
|
|
PRECONDITION(mol, "no molecule");
|
|
|
|
ROMol::AtomIterator atomIt;
|
|
for (atomIt = mol->beginAtoms(); atomIt != mol->endAtoms(); atomIt++) {
|
|
transformAtom(*atomIt, tform);
|
|
}
|
|
}
|
|
|
|
RDGeom::Point3D computeCentroid(const Conformer &conf, bool ignoreHs,
|
|
const std::vector<double> *weights) {
|
|
PRECONDITION(!weights || weights->size() >= conf.getNumAtoms(),
|
|
"bad weights vector");
|
|
RDGeom::Point3D res(0.0, 0.0, 0.0);
|
|
const ROMol &mol = conf.getOwningMol();
|
|
double wSum = 0.0;
|
|
for (unsigned int i = 0; i < conf.getNumAtoms(); ++i) {
|
|
if (ignoreHs && mol.getAtomWithIdx(i)->getAtomicNum() == 1) {
|
|
continue;
|
|
}
|
|
double w = (weights ? weights->at(i) : 1.0);
|
|
wSum += w;
|
|
res += conf.getAtomPos(i) * w;
|
|
}
|
|
res /= wSum;
|
|
return res;
|
|
}
|
|
|
|
namespace {
|
|
void computeCovarianceTerms(const Conformer &conf,
|
|
const RDGeom::Point3D ¢er, double &xx,
|
|
double &xy, double &xz, double &yy, double &yz,
|
|
double &zz, bool normalize, bool ignoreHs,
|
|
const std::vector<double> *weights) {
|
|
PRECONDITION(!weights || weights->size() >= conf.getNumAtoms(),
|
|
"bad weights vector");
|
|
|
|
xx = xy = xz = yy = yz = zz = 0.0;
|
|
const ROMol &mol = conf.getOwningMol();
|
|
double wSum = 0.0;
|
|
for (ROMol::ConstAtomIterator cai = mol.beginAtoms(); cai != mol.endAtoms();
|
|
cai++) {
|
|
if (((*cai)->getAtomicNum() == 1) && (ignoreHs)) {
|
|
continue;
|
|
}
|
|
RDGeom::Point3D loc = conf.getAtomPos((*cai)->getIdx());
|
|
loc -= center;
|
|
double w = 1.0;
|
|
if (weights) {
|
|
w = (*weights)[(*cai)->getIdx()];
|
|
}
|
|
wSum += w;
|
|
xx += w * loc.x * loc.x;
|
|
xy += w * loc.x * loc.y;
|
|
xz += w * loc.x * loc.z;
|
|
yy += w * loc.y * loc.y;
|
|
yz += w * loc.y * loc.z;
|
|
zz += w * loc.z * loc.z;
|
|
}
|
|
if (normalize) {
|
|
xx /= wSum;
|
|
xy /= wSum;
|
|
xz /= wSum;
|
|
yy /= wSum;
|
|
yz /= wSum;
|
|
zz /= wSum;
|
|
}
|
|
}
|
|
|
|
#ifndef RDK_HAS_EIGEN3
|
|
RDNumeric::DoubleSymmMatrix *computeCovarianceMatrix(
|
|
const Conformer &conf, const RDGeom::Point3D ¢er, bool normalize,
|
|
bool ignoreHs) {
|
|
double xx, xy, xz, yy, yz, zz;
|
|
computeCovarianceTerms(conf, center, xx, xy, xz, yy, yz, zz, normalize,
|
|
ignoreHs, nullptr);
|
|
auto *res = new RDNumeric::DoubleSymmMatrix(3, 3);
|
|
res->setVal(0, 0, xx);
|
|
res->setVal(0, 1, xy);
|
|
res->setVal(0, 2, xz);
|
|
res->setVal(1, 1, yy);
|
|
res->setVal(1, 2, yz);
|
|
res->setVal(2, 2, zz);
|
|
return res;
|
|
}
|
|
#endif
|
|
|
|
void computeInertiaTerms(const Conformer &conf, const RDGeom::Point3D ¢er,
|
|
double &xx, double &xy, double &xz, double &yy,
|
|
double &yz, double &zz, bool ignoreHs,
|
|
const std::vector<double> *weights) {
|
|
PRECONDITION(!weights || weights->size() >= conf.getNumAtoms(),
|
|
"bad weights vector");
|
|
|
|
xx = xy = xz = yy = yz = zz = 0.0;
|
|
const ROMol &mol = conf.getOwningMol();
|
|
for (ROMol::ConstAtomIterator cai = mol.beginAtoms(); cai != mol.endAtoms();
|
|
cai++) {
|
|
if (((*cai)->getAtomicNum() == 1) && (ignoreHs)) {
|
|
continue;
|
|
}
|
|
RDGeom::Point3D loc = conf.getAtomPos((*cai)->getIdx());
|
|
loc -= center;
|
|
double w = 1.0;
|
|
if (weights) {
|
|
w = (*weights)[(*cai)->getIdx()];
|
|
}
|
|
xx += w * (loc.y * loc.y + loc.z * loc.z);
|
|
yy += w * (loc.x * loc.x + loc.z * loc.z);
|
|
zz += w * (loc.y * loc.y + loc.x * loc.x);
|
|
xy -= w * loc.x * loc.y;
|
|
xz -= w * loc.x * loc.z;
|
|
yz -= w * loc.z * loc.y;
|
|
}
|
|
}
|
|
} // namespace
|
|
|
|
#ifdef RDK_HAS_EIGEN3
|
|
#include <Eigen/Dense>
|
|
|
|
namespace {
|
|
bool getEigenValEigenVectHelper(Eigen::Matrix3d &eigVecs,
|
|
Eigen::Vector3d &eigVals, double sumXX,
|
|
double sumXY, double sumXZ, double sumYY,
|
|
double sumYZ, double sumZZ) {
|
|
Eigen::Matrix3d mat;
|
|
mat << sumXX, sumXY, sumXZ, sumXY, sumYY, sumYZ, sumXZ, sumYZ, sumZZ;
|
|
// std::cerr<<"getEigenValEigenVectHelper matrix: "<<mat<<std::endl;
|
|
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(mat);
|
|
if (eigensolver.info() != Eigen::Success) {
|
|
BOOST_LOG(rdErrorLog) << "eigenvalue calculation did not converge"
|
|
<< std::endl;
|
|
return false;
|
|
}
|
|
eigVecs = eigensolver.eigenvectors();
|
|
// std::cerr<<"getEigenValEigenVectHelper eigVecs: "<<eigVecs<<std::endl;
|
|
eigVals = eigensolver.eigenvalues();
|
|
// std::cerr<<"getEigenValEigenVectHelper eigVals: "<<eigVals<<std::endl;
|
|
return true;
|
|
}
|
|
|
|
bool getEigenValEigenVectFromCovMat(const RDKit::Conformer &conf,
|
|
Eigen::Matrix3d &eigVecs,
|
|
Eigen::Vector3d &eigVals,
|
|
const RDGeom::Point3D &origin,
|
|
bool ignoreHs, bool normalizeCovar,
|
|
const std::vector<double> *weights) {
|
|
PRECONDITION((!weights || weights->size() >= conf.getNumAtoms()),
|
|
"bad weights vector");
|
|
|
|
// std::cerr << "getEigenValEigenVectFromCovMat ignoreHs " << ignoreHs << "
|
|
// normalizeCovar " << normalizeCovar << " weights " << weights << " origin "
|
|
// << origin.x << "," << origin.y << "," << origin.z << std::endl;
|
|
double sumXX, sumXY, sumXZ, sumYY, sumYZ, sumZZ;
|
|
computeCovarianceTerms(conf, origin, sumXX, sumXY, sumXZ, sumYY, sumYZ, sumZZ,
|
|
normalizeCovar, ignoreHs, weights);
|
|
|
|
return getEigenValEigenVectHelper(eigVecs, eigVals, sumXX, sumXY, sumXZ,
|
|
sumYY, sumYZ, sumZZ);
|
|
}
|
|
} // namespace
|
|
|
|
bool computePrincipalAxesAndMoments(const RDKit::Conformer &conf,
|
|
Eigen::Matrix3d &axes,
|
|
Eigen::Vector3d &moments, bool ignoreHs,
|
|
bool force,
|
|
const std::vector<double> *weights) {
|
|
PRECONDITION((!weights || weights->size() >= conf.getNumAtoms()),
|
|
"bad weights vector");
|
|
const char *axesPropName = ignoreHs ? "_principalAxes_noH" : "_principalAxes";
|
|
const char *momentsPropName =
|
|
ignoreHs ? "_principalMoments_noH" : "_principalMoments";
|
|
const ROMol &mol = conf.getOwningMol();
|
|
if (!weights && !force && mol.hasProp(axesPropName) &&
|
|
mol.hasProp(momentsPropName)) {
|
|
mol.getProp(axesPropName, axes);
|
|
mol.getProp(momentsPropName, moments);
|
|
return true;
|
|
}
|
|
auto origin = computeCentroid(conf, ignoreHs, weights);
|
|
|
|
double sumXX, sumXY, sumXZ, sumYY, sumYZ, sumZZ;
|
|
computeInertiaTerms(conf, origin, sumXX, sumXY, sumXZ, sumYY, sumYZ, sumZZ,
|
|
ignoreHs, weights);
|
|
|
|
if (!getEigenValEigenVectHelper(axes, moments, sumXX, sumXY, sumXZ, sumYY,
|
|
sumYZ, sumZZ)) {
|
|
return false;
|
|
}
|
|
|
|
if (!weights) {
|
|
mol.setProp(axesPropName, axes, true);
|
|
mol.setProp(momentsPropName, moments, true);
|
|
}
|
|
return true;
|
|
}
|
|
|
|
bool computePrincipalAxesAndMomentsFromGyrationMatrix(
|
|
const RDKit::Conformer &conf, Eigen::Matrix3d &axes,
|
|
Eigen::Vector3d &moments, bool ignoreHs, bool force,
|
|
const std::vector<double> *weights) {
|
|
const char *axesPropName =
|
|
ignoreHs ? "_principalAxes_noH_cov" : "_principalAxes_cov";
|
|
const char *momentsPropName =
|
|
ignoreHs ? "_principalMoments_noH_cov" : "_principalMoments_cov";
|
|
const ROMol &mol = conf.getOwningMol();
|
|
if (!weights && !force && mol.hasProp(axesPropName) &&
|
|
mol.hasProp(momentsPropName)) {
|
|
mol.getProp(axesPropName, axes);
|
|
mol.getProp(momentsPropName, moments);
|
|
return true;
|
|
}
|
|
auto origin = computeCentroid(conf, ignoreHs, weights);
|
|
bool res = getEigenValEigenVectFromCovMat(conf, axes, moments, origin,
|
|
ignoreHs, true, weights);
|
|
if (res && !weights) {
|
|
conf.getOwningMol().setProp(axesPropName, axes, true);
|
|
conf.getOwningMol().setProp(momentsPropName, moments, true);
|
|
}
|
|
return res;
|
|
}
|
|
|
|
RDGeom::Transform3D *computeCanonicalTransform(const Conformer &conf,
|
|
const RDGeom::Point3D *center,
|
|
bool normalizeCovar,
|
|
bool ignoreHs) {
|
|
constexpr unsigned int DIM = 3;
|
|
RDGeom::Point3D origin;
|
|
if (!center) {
|
|
origin = computeCentroid(conf, ignoreHs);
|
|
} else {
|
|
origin = (*center);
|
|
}
|
|
unsigned int nAtms = conf.getNumAtoms();
|
|
auto *trans = new RDGeom::Transform3D;
|
|
trans->setToIdentity();
|
|
|
|
// if we have a single atom system we don't need to do anyhting setting
|
|
// translation is sufficient
|
|
if (nAtms > 1) {
|
|
Eigen::Matrix3d eigVecs;
|
|
Eigen::Vector3d eigVals;
|
|
if (getEigenValEigenVectFromCovMat(conf, eigVecs, eigVals, origin, ignoreHs,
|
|
normalizeCovar, nullptr)) {
|
|
std::vector<std::pair<unsigned int, double>> eigValsSorted;
|
|
CHECK_INVARIANT(eigVals.size() == DIM, "less eigenvalues than expected");
|
|
eigValsSorted.reserve(DIM);
|
|
for (unsigned int i = 0; i < DIM; ++i) {
|
|
eigValsSorted.emplace_back(i, eigVals(i));
|
|
}
|
|
std::sort(eigValsSorted.begin(), eigValsSorted.end(),
|
|
[](const std::pair<unsigned int, double> &a,
|
|
const std::pair<unsigned int, double> &b) {
|
|
return (a.second > b.second);
|
|
});
|
|
for (unsigned int col = 0; col < DIM; ++col) {
|
|
unsigned int colSorted = eigValsSorted.at(col).first;
|
|
for (unsigned int row = 0; row < DIM; ++row) {
|
|
trans->setVal(col, row, eigVecs(row, colSorted));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
origin *= -1.0;
|
|
|
|
// Check that we have received a right-handed coordinate axis -
|
|
// this is not guaranteed from the jacobi routine. The test
|
|
// is to see whether the determinant of the eigenvector matrix (which
|
|
// is the equivalent of ( {vec1 cross vec2} dot vec3 ) is +1 or
|
|
// -1. In the latter case, it is left-handed, so reverse the sign
|
|
// of one of the vectors (in this case the z). This fix came from
|
|
// Andy Grant, RIP.
|
|
double test =
|
|
trans->getVal(0, 2) * (trans->getVal(1, 0) * trans->getVal(2, 1) -
|
|
trans->getVal(2, 0) * trans->getVal(1, 1)) +
|
|
trans->getVal(1, 2) * (trans->getVal(2, 0) * trans->getVal(0, 1) -
|
|
trans->getVal(2, 1) * trans->getVal(0, 0)) +
|
|
trans->getVal(2, 2) * (trans->getVal(0, 0) * trans->getVal(1, 1) -
|
|
trans->getVal(1, 0) * trans->getVal(0, 1));
|
|
if (test < 0.0) {
|
|
for (auto i = 0; i < 3; ++i) {
|
|
trans->setVal(2, i, trans->getVal(2, i) * -1);
|
|
}
|
|
}
|
|
trans->TransformPoint(origin);
|
|
trans->SetTranslation(origin);
|
|
|
|
return trans;
|
|
}
|
|
#else
|
|
RDGeom::Transform3D *computeCanonicalTransform(const Conformer &conf,
|
|
const RDGeom::Point3D *center,
|
|
bool normalizeCovar,
|
|
bool ignoreHs) {
|
|
RDGeom::Point3D origin;
|
|
if (!center) {
|
|
origin = computeCentroid(conf, ignoreHs);
|
|
} else {
|
|
origin = (*center);
|
|
}
|
|
RDNumeric::DoubleSymmMatrix *covMat =
|
|
computeCovarianceMatrix(conf, origin, normalizeCovar, ignoreHs);
|
|
// find the eigen values and eigen vectors for the covMat
|
|
RDNumeric::DoubleMatrix eigVecs(3, 3);
|
|
RDNumeric::DoubleVector eigVals(3);
|
|
// if we have a single atom system we don't need to do anyhting other than
|
|
// setting translation
|
|
// translation
|
|
unsigned int nAtms = conf.getNumAtoms();
|
|
auto *trans = new RDGeom::Transform3D;
|
|
|
|
// set the translation
|
|
origin *= -1.0;
|
|
// trans->SetTranslation(origin);
|
|
// if we have a single atom system we don't need to do anyhting setting
|
|
// translation is sufficient
|
|
if (nAtms > 1) {
|
|
RDNumeric::EigenSolvers::powerEigenSolver(3, *covMat, eigVals, eigVecs,
|
|
conf.getNumAtoms());
|
|
// deal with zero eigen value systems
|
|
unsigned int i, j, dim = 3;
|
|
for (i = 0; i < 3; ++i) {
|
|
if (fabs(eigVals.getVal(i)) < EIGEN_TOLERANCE) {
|
|
dim--;
|
|
}
|
|
}
|
|
CHECK_INVARIANT(dim >= 1, "");
|
|
if (dim < 3) {
|
|
RDGeom::Point3D first(eigVecs.getVal(0, 0), eigVecs.getVal(0, 1),
|
|
eigVecs.getVal(0, 2));
|
|
if (dim == 1) {
|
|
// pick an arbitrary eigen vector perpendicular to the first vector
|
|
RDGeom::Point3D second(first.getPerpendicular());
|
|
eigVecs.setVal(1, 0, second.x);
|
|
eigVecs.setVal(1, 1, second.y);
|
|
eigVecs.setVal(1, 2, second.z);
|
|
if (eigVals.getVal(0) > 1.0) {
|
|
eigVals.setVal(1, 1.0);
|
|
} else {
|
|
eigVals.setVal(1, eigVals.getVal(0) / 2.0);
|
|
}
|
|
}
|
|
RDGeom::Point3D second(eigVecs.getVal(1, 0), eigVecs.getVal(1, 1),
|
|
eigVecs.getVal(1, 2));
|
|
// pick the third eigen vector perpendicular to the first two
|
|
RDGeom::Point3D third = first.crossProduct(second);
|
|
eigVecs.setVal(2, 0, third.x);
|
|
eigVecs.setVal(2, 1, third.y);
|
|
eigVecs.setVal(2, 2, third.z);
|
|
if (eigVals.getVal(1) > 1.0) {
|
|
eigVals.setVal(2, 1.0);
|
|
} else {
|
|
eigVals.setVal(2, eigVals.getVal(1) / 2.0);
|
|
}
|
|
}
|
|
// now set the transformation
|
|
for (i = 0; i < 3; ++i) {
|
|
for (j = 0; j < 3; ++j) {
|
|
trans->setVal(i, j, eigVecs.getVal(i, j));
|
|
}
|
|
}
|
|
} // end of multiple atom system
|
|
|
|
// In some situations we can end up with one or more negative values on the
|
|
// diagonal. An odd number of these will result in an inversion of the
|
|
// structure, so we need to check for that and, if necessary, correct by
|
|
// negating one row
|
|
if (trans->getVal(0, 0) * trans->getVal(1, 1) * trans->getVal(2, 2) < 0) {
|
|
for (auto i = 0; i < 3; ++i) {
|
|
trans->setVal(2, i, trans->getVal(2, i) * -1);
|
|
}
|
|
}
|
|
trans->TransformPoint(origin);
|
|
|
|
trans->SetTranslation(origin);
|
|
delete covMat;
|
|
|
|
return trans;
|
|
}
|
|
#endif
|
|
|
|
void canonicalizeConformer(Conformer &conf, const RDGeom::Point3D *center,
|
|
bool normalizeCovar, bool ignoreHs) {
|
|
RDGeom::Transform3D *trans =
|
|
computeCanonicalTransform(conf, center, normalizeCovar, ignoreHs);
|
|
transformConformer(conf, *trans);
|
|
delete trans;
|
|
}
|
|
|
|
void canonicalizeMol(RDKit::ROMol &mol, bool normalizeCovar, bool ignoreHs) {
|
|
ROMol::ConformerIterator ci;
|
|
for (ci = mol.beginConformers(); ci != mol.endConformers(); ci++) {
|
|
canonicalizeConformer(*(*ci), nullptr, normalizeCovar, ignoreHs);
|
|
}
|
|
}
|
|
|
|
void transformConformer(Conformer &conf, const RDGeom::Transform3D &trans) {
|
|
RDGeom::POINT3D_VECT &positions = conf.getPositions();
|
|
RDGeom::POINT3D_VECT_I pi;
|
|
for (pi = positions.begin(); pi != positions.end(); ++pi) {
|
|
trans.TransformPoint(*pi);
|
|
}
|
|
}
|
|
|
|
void transformMolSubstanceGroups(ROMol &mol, const RDGeom::Transform3D &trans) {
|
|
auto &sgs = getSubstanceGroups(mol);
|
|
for (auto &sg : sgs) {
|
|
for (auto &brk : sg.getBrackets()) {
|
|
trans.TransformPoint(brk[0]);
|
|
trans.TransformPoint(brk[1]);
|
|
trans.TransformPoint(brk[2]);
|
|
}
|
|
for (auto &cs : sg.getCStates()) {
|
|
trans.TransformPoint(cs.vector);
|
|
}
|
|
}
|
|
}
|
|
namespace {
|
|
void _toBeMovedIdxList(const ROMol &mol, unsigned int iAtomId,
|
|
unsigned int jAtomId, std::list<unsigned int> &alist) {
|
|
unsigned int nAtoms = mol.getNumAtoms();
|
|
boost::dynamic_bitset<> visitedIdx(nAtoms);
|
|
std::stack<unsigned int> stack;
|
|
stack.push(jAtomId);
|
|
visitedIdx[iAtomId] = 1;
|
|
visitedIdx[jAtomId] = 1;
|
|
unsigned int tIdx;
|
|
unsigned int wIdx;
|
|
ROMol::ADJ_ITER nbrIdx;
|
|
ROMol::ADJ_ITER endNbrs;
|
|
bool doMainLoop;
|
|
while (stack.size()) {
|
|
doMainLoop = false;
|
|
tIdx = stack.top();
|
|
const Atom *tAtom = mol.getAtomWithIdx(tIdx);
|
|
boost::tie(nbrIdx, endNbrs) = mol.getAtomNeighbors(tAtom);
|
|
unsigned int eIdx;
|
|
for (eIdx = 0; nbrIdx != endNbrs; ++nbrIdx, ++eIdx) {
|
|
wIdx = (mol[*nbrIdx])->getIdx();
|
|
if (!visitedIdx[wIdx]) {
|
|
visitedIdx[wIdx] = 1;
|
|
stack.push(wIdx);
|
|
doMainLoop = true;
|
|
break;
|
|
}
|
|
}
|
|
if (doMainLoop) {
|
|
continue;
|
|
}
|
|
visitedIdx[tIdx] = 1;
|
|
stack.pop();
|
|
}
|
|
alist.clear();
|
|
for (unsigned int i = 0; i < nAtoms; ++i) {
|
|
if (visitedIdx[i] && (i != iAtomId)) {
|
|
alist.push_back(i);
|
|
}
|
|
}
|
|
}
|
|
} // namespace
|
|
double getBondLength(const Conformer &conf, unsigned int iAtomId,
|
|
unsigned int jAtomId) {
|
|
const RDGeom::POINT3D_VECT &pos = conf.getPositions();
|
|
URANGE_CHECK(iAtomId, pos.size());
|
|
URANGE_CHECK(jAtomId, pos.size());
|
|
|
|
return (pos[iAtomId] - pos[jAtomId]).length();
|
|
}
|
|
|
|
void setBondLength(Conformer &conf, unsigned int iAtomId, unsigned int jAtomId,
|
|
double value) {
|
|
RDGeom::POINT3D_VECT &pos = conf.getPositions();
|
|
URANGE_CHECK(iAtomId, pos.size());
|
|
URANGE_CHECK(jAtomId, pos.size());
|
|
ROMol &mol = conf.getOwningMol();
|
|
Bond *bond = mol.getBondBetweenAtoms(iAtomId, jAtomId);
|
|
if (!bond) {
|
|
throw ValueErrorException("atoms i and j must be bonded");
|
|
}
|
|
if (queryIsBondInRing(bond)) {
|
|
throw ValueErrorException("bond (i,j) must not belong to a ring");
|
|
}
|
|
RDGeom::Point3D v = pos[iAtomId] - pos[jAtomId];
|
|
double origValue = v.length();
|
|
if (origValue <= 1.e-8) {
|
|
throw ValueErrorException("atoms i and j have identical 3D coordinates");
|
|
}
|
|
|
|
// get all atoms bonded to j
|
|
std::list<unsigned int> alist;
|
|
_toBeMovedIdxList(mol, iAtomId, jAtomId, alist);
|
|
v *= (value / origValue - 1.);
|
|
for (unsigned int &it : alist) {
|
|
pos[it] -= v;
|
|
}
|
|
}
|
|
|
|
double getAngleRad(const Conformer &conf, unsigned int iAtomId,
|
|
unsigned int jAtomId, unsigned int kAtomId) {
|
|
const RDGeom::POINT3D_VECT &pos = conf.getPositions();
|
|
URANGE_CHECK(iAtomId, pos.size());
|
|
URANGE_CHECK(jAtomId, pos.size());
|
|
URANGE_CHECK(kAtomId, pos.size());
|
|
RDGeom::Point3D rJI = pos[iAtomId] - pos[jAtomId];
|
|
double rJISqLength = rJI.lengthSq();
|
|
if (rJISqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms i and j have identical 3D coordinates");
|
|
}
|
|
RDGeom::Point3D rJK = pos[kAtomId] - pos[jAtomId];
|
|
double rJKSqLength = rJK.lengthSq();
|
|
if (rJKSqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms j and k have identical 3D coordinates");
|
|
}
|
|
return rJI.angleTo(rJK);
|
|
}
|
|
|
|
void setAngleRad(Conformer &conf, unsigned int iAtomId, unsigned int jAtomId,
|
|
unsigned int kAtomId, double value) {
|
|
RDGeom::POINT3D_VECT &pos = conf.getPositions();
|
|
URANGE_CHECK(iAtomId, pos.size());
|
|
URANGE_CHECK(jAtomId, pos.size());
|
|
URANGE_CHECK(kAtomId, pos.size());
|
|
ROMol &mol = conf.getOwningMol();
|
|
Bond *bondJI = mol.getBondBetweenAtoms(jAtomId, iAtomId);
|
|
if (!bondJI) {
|
|
throw ValueErrorException("atoms i and j must be bonded");
|
|
}
|
|
Bond *bondJK = mol.getBondBetweenAtoms(jAtomId, kAtomId);
|
|
if (!bondJK) {
|
|
throw ValueErrorException("atoms j and k must be bonded");
|
|
}
|
|
if (queryIsBondInRing(bondJI) && queryIsBondInRing(bondJK)) {
|
|
throw ValueErrorException(
|
|
"bonds (i,j) and (j,k) must not both belong to a ring");
|
|
}
|
|
|
|
RDGeom::Point3D rJI = pos[iAtomId] - pos[jAtomId];
|
|
double rJISqLength = rJI.lengthSq();
|
|
if (rJISqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms i and j have identical 3D coordinates");
|
|
}
|
|
RDGeom::Point3D rJK = pos[kAtomId] - pos[jAtomId];
|
|
double rJKSqLength = rJK.lengthSq();
|
|
if (rJKSqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms j and k have identical 3D coordinates");
|
|
}
|
|
|
|
// we only need to rotate by delta with respect to the current angle value
|
|
value -= rJI.angleTo(rJK);
|
|
RDGeom::Point3D &rotAxisBegin = pos[jAtomId];
|
|
// our rotation axis is the normal to the plane of atoms i, j, k
|
|
RDGeom::Point3D rotAxisEnd = rJI.crossProduct(rJK) + pos[jAtomId];
|
|
RDGeom::Point3D rotAxis = rotAxisEnd - rotAxisBegin;
|
|
rotAxis.normalize();
|
|
// get all atoms bonded to j and loop through them
|
|
std::list<unsigned int> alist;
|
|
_toBeMovedIdxList(mol, jAtomId, kAtomId, alist);
|
|
for (unsigned int &it : alist) {
|
|
// translate atom so that it coincides with the origin of rotation
|
|
pos[it] -= rotAxisBegin;
|
|
// rotate around our rotation axis
|
|
RDGeom::Transform3D rotByAngle;
|
|
rotByAngle.SetRotation(value, rotAxis);
|
|
rotByAngle.TransformPoint(pos[it]);
|
|
// translate atom back
|
|
pos[it] += rotAxisBegin;
|
|
}
|
|
}
|
|
|
|
double getDihedralRad(const Conformer &conf, unsigned int iAtomId,
|
|
unsigned int jAtomId, unsigned int kAtomId,
|
|
unsigned int lAtomId) {
|
|
const RDGeom::POINT3D_VECT &pos = conf.getPositions();
|
|
URANGE_CHECK(iAtomId, pos.size());
|
|
URANGE_CHECK(jAtomId, pos.size());
|
|
URANGE_CHECK(kAtomId, pos.size());
|
|
URANGE_CHECK(lAtomId, pos.size());
|
|
RDGeom::Point3D rIJ = pos[jAtomId] - pos[iAtomId];
|
|
double rIJSqLength = rIJ.lengthSq();
|
|
if (rIJSqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms i and j have identical 3D coordinates");
|
|
}
|
|
RDGeom::Point3D rJK = pos[kAtomId] - pos[jAtomId];
|
|
double rJKSqLength = rJK.lengthSq();
|
|
if (rJKSqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms j and k have identical 3D coordinates");
|
|
}
|
|
RDGeom::Point3D rKL = pos[lAtomId] - pos[kAtomId];
|
|
double rKLSqLength = rKL.lengthSq();
|
|
if (rKLSqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms k and l have identical 3D coordinates");
|
|
}
|
|
|
|
RDGeom::Point3D nIJK = rIJ.crossProduct(rJK);
|
|
double nIJKSqLength = nIJK.lengthSq();
|
|
RDGeom::Point3D nJKL = rJK.crossProduct(rKL);
|
|
double nJKLSqLength = nJKL.lengthSq();
|
|
RDGeom::Point3D m = nIJK.crossProduct(rJK);
|
|
// we want a signed dihedral, that's why we use atan2 instead of acos
|
|
return -atan2(m.dotProduct(nJKL) / sqrt(nJKLSqLength * m.lengthSq()),
|
|
nIJK.dotProduct(nJKL) / sqrt(nIJKSqLength * nJKLSqLength));
|
|
}
|
|
|
|
void setDihedralRad(Conformer &conf, unsigned int iAtomId, unsigned int jAtomId,
|
|
unsigned int kAtomId, unsigned int lAtomId, double value) {
|
|
RDGeom::POINT3D_VECT &pos = conf.getPositions();
|
|
URANGE_CHECK(iAtomId, pos.size());
|
|
URANGE_CHECK(jAtomId, pos.size());
|
|
URANGE_CHECK(kAtomId, pos.size());
|
|
URANGE_CHECK(lAtomId, pos.size());
|
|
ROMol &mol = conf.getOwningMol();
|
|
Bond *bondJK = mol.getBondBetweenAtoms(jAtomId, kAtomId);
|
|
if (!bondJK) {
|
|
throw ValueErrorException("atoms j and k must be bonded");
|
|
}
|
|
|
|
if (queryIsBondInRing(bondJK)) {
|
|
throw ValueErrorException("bond (j,k) must not belong to a ring");
|
|
}
|
|
RDGeom::Point3D rIJ = pos[jAtomId] - pos[iAtomId];
|
|
double rIJSqLength = rIJ.lengthSq();
|
|
if (rIJSqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms i and j have identical 3D coordinates");
|
|
}
|
|
RDGeom::Point3D rJK = pos[kAtomId] - pos[jAtomId];
|
|
double rJKSqLength = rJK.lengthSq();
|
|
if (rJKSqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms j and k have identical 3D coordinates");
|
|
}
|
|
RDGeom::Point3D rKL = pos[lAtomId] - pos[kAtomId];
|
|
double rKLSqLength = rKL.lengthSq();
|
|
if (rKLSqLength <= 1.e-16) {
|
|
throw ValueErrorException("atoms k and l have identical 3D coordinates");
|
|
}
|
|
|
|
RDGeom::Point3D nIJK = rIJ.crossProduct(rJK);
|
|
double nIJKSqLength = nIJK.lengthSq();
|
|
RDGeom::Point3D nJKL = rJK.crossProduct(rKL);
|
|
double nJKLSqLength = nJKL.lengthSq();
|
|
RDGeom::Point3D m = nIJK.crossProduct(rJK);
|
|
// we only need to rotate by delta with respect to the current dihedral value
|
|
value -= -atan2(m.dotProduct(nJKL) / sqrt(nJKLSqLength * m.lengthSq()),
|
|
nIJK.dotProduct(nJKL) / sqrt(nIJKSqLength * nJKLSqLength));
|
|
// our rotation axis is the (j,k) bond
|
|
RDGeom::Point3D &rotAxisBegin = pos[jAtomId];
|
|
RDGeom::Point3D &rotAxisEnd = pos[kAtomId];
|
|
RDGeom::Point3D rotAxis = rotAxisEnd - rotAxisBegin;
|
|
rotAxis.normalize();
|
|
// get all atoms bonded to k and loop through them
|
|
std::list<unsigned int> alist;
|
|
_toBeMovedIdxList(mol, jAtomId, kAtomId, alist);
|
|
for (unsigned int &it : alist) {
|
|
// translate atom so that it coincides with the origin of rotation
|
|
pos[it] -= rotAxisBegin;
|
|
// rotate around our rotation axis
|
|
RDGeom::Transform3D rotByAngle;
|
|
rotByAngle.SetRotation(value, rotAxis);
|
|
rotByAngle.TransformPoint(pos[it]);
|
|
// translate atom back
|
|
pos[it] += rotAxisBegin;
|
|
}
|
|
}
|
|
} // namespace MolTransforms
|