mirror of
https://github.com/rdkit/rdkit.git
synced 2026-06-04 21:54:27 +08:00
756 lines
28 KiB
C++
756 lines
28 KiB
C++
// $Id$
|
|
//
|
|
// Copyright (C) 2004-2012 Greg Landrum and Rational Discovery LLC
|
|
//
|
|
// @@ 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 "Embedder.h"
|
|
#include <DistGeom/BoundsMatrix.h>
|
|
#include <DistGeom/DistGeomUtils.h>
|
|
#include <DistGeom/TriangleSmooth.h>
|
|
#include <DistGeom/ChiralViolationContrib.h>
|
|
#include "BoundsMatrixBuilder.h"
|
|
#include <ForceField/ForceField.h>
|
|
#include <GraphMol/ROMol.h>
|
|
#include <GraphMol/Atom.h>
|
|
#include <GraphMol/AtomIterators.h>
|
|
|
|
#include <GraphMol/Conformer.h>
|
|
#include <RDGeneral/types.h>
|
|
#include <RDGeneral/RDLog.h>
|
|
#include <RDGeneral/Exceptions.h>
|
|
|
|
#include <Geometry/Transform3D.h>
|
|
#include <Numerics/Alignment/AlignPoints.h>
|
|
#include <DistGeom/ChiralSet.h>
|
|
#include <GraphMol/MolOps.h>
|
|
#include <GraphMol/ForceFieldHelpers/CrystalFF/TorsionPreferences.h>
|
|
#include <boost/dynamic_bitset.hpp>
|
|
#include <iomanip>
|
|
#include <RDGeneral/RDThreads.h>
|
|
|
|
#define ERROR_TOL 0.00001
|
|
|
|
namespace RDKit {
|
|
namespace DGeomHelpers {
|
|
typedef std::pair<int, int> INT_PAIR;
|
|
typedef std::vector<INT_PAIR> INT_PAIR_VECT;
|
|
|
|
bool _sameSide(const RDGeom::Point3D &v1, const RDGeom::Point3D &v2,
|
|
const RDGeom::Point3D &v3, const RDGeom::Point3D &v4,
|
|
const RDGeom::Point3D &p0, double tol = 0.1) {
|
|
RDGeom::Point3D normal = (v2 - v1).crossProduct(v3 - v1);
|
|
double d1 = normal.dotProduct(v4 - v1);
|
|
double d2 = normal.dotProduct(p0 - v1);
|
|
// std::cerr<<" "<<d1<<" - " <<d2<<std::endl;
|
|
if (fabs(d1) < tol || fabs(d2) < tol) return false;
|
|
return !((d1 < 0.) ^ (d2 < 0.));
|
|
}
|
|
bool _centerInVolume(const DistGeom::ChiralSetPtr &chiralSet,
|
|
const RDGeom::PointPtrVect &positions) {
|
|
if (chiralSet->d_idx0 ==
|
|
chiralSet->d_idx4) { // this happens for three-coordinate centers
|
|
return true;
|
|
}
|
|
|
|
RDGeom::Point3D p0((*positions[chiralSet->d_idx0])[0],
|
|
(*positions[chiralSet->d_idx0])[1],
|
|
(*positions[chiralSet->d_idx0])[2]);
|
|
RDGeom::Point3D p1((*positions[chiralSet->d_idx1])[0],
|
|
(*positions[chiralSet->d_idx1])[1],
|
|
(*positions[chiralSet->d_idx1])[2]);
|
|
RDGeom::Point3D p2((*positions[chiralSet->d_idx2])[0],
|
|
(*positions[chiralSet->d_idx2])[1],
|
|
(*positions[chiralSet->d_idx2])[2]);
|
|
RDGeom::Point3D p3((*positions[chiralSet->d_idx3])[0],
|
|
(*positions[chiralSet->d_idx3])[1],
|
|
(*positions[chiralSet->d_idx3])[2]);
|
|
RDGeom::Point3D p4((*positions[chiralSet->d_idx4])[0],
|
|
(*positions[chiralSet->d_idx4])[1],
|
|
(*positions[chiralSet->d_idx4])[2]);
|
|
// RDGeom::Point3D centroid = (p1+p2+p3+p4)/4.;
|
|
|
|
bool res = _sameSide(p1, p2, p3, p4, p0) && _sameSide(p2, p3, p4, p1, p0) &&
|
|
_sameSide(p3, p4, p1, p2, p0) && _sameSide(p4, p1, p2, p3, p0);
|
|
// std::cerr<<"civ:"<<chiralSet->d_idx0<<" "<<chiralSet->d_idx1<<"
|
|
// "<<chiralSet->d_idx2<<" "<<chiralSet->d_idx3<<"
|
|
// "<<chiralSet->d_idx4<<"->"<<res<<"|"<<std::endl;
|
|
return res;
|
|
}
|
|
bool _boundsFulfilled(const std::vector<int> &atoms,
|
|
const DistGeom::BoundsMatrix &mmat,
|
|
const RDGeom::PointPtrVect &positions) {
|
|
// unsigned int N = mmat.numRows();
|
|
// std::cerr << N << " " << atoms.size() << std::endl;
|
|
// loop over all pair of atoms
|
|
for (unsigned int i = 0; i < atoms.size() - 1; ++i) {
|
|
for (unsigned int j = i + 1; j < atoms.size(); ++j) {
|
|
int a1 = atoms[i];
|
|
int a2 = atoms[j];
|
|
RDGeom::Point3D p0((*positions[a1])[0], (*positions[a1])[1],
|
|
(*positions[a1])[2]);
|
|
RDGeom::Point3D p1((*positions[a2])[0], (*positions[a2])[1],
|
|
(*positions[a2])[2]);
|
|
double d2 = (p0 - p1).length(); // distance
|
|
double lb = mmat.getLowerBound(a1, a2);
|
|
double ub = mmat.getUpperBound(a1, a2); // bounds
|
|
if (((d2 < lb) && (fabs(d2 - lb) > 0.17)) ||
|
|
((d2 > ub) && (fabs(d2 - ub) > 0.17))) {
|
|
// std::cerr << a1 << " " << a2 << ":" << d2 << " " << lb << " " << ub
|
|
// << " " << fabs(d2-lb) << " " << fabs(d2-ub) << std::endl;
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
// the minimization using experimental torsion angle preferences
|
|
void _minimizeWithExpTorsions(
|
|
RDGeom::PointPtrVect &positions, DistGeom::BoundsMatPtr mmat,
|
|
double optimizerForceTol, double basinThresh,
|
|
const std::vector<std::pair<int, int> > &bonds,
|
|
const std::vector<std::vector<int> > &angles,
|
|
const std::vector<std::vector<int> > &expTorsionAtoms,
|
|
const std::vector<std::pair<std::vector<int>, std::vector<double> > >
|
|
&expTorsionAngles,
|
|
const std::vector<std::vector<int> > &improperAtoms,
|
|
const std::vector<int> &atomNums, bool useBasicKnowledge) {
|
|
RDUNUSED_PARAM(basinThresh);
|
|
|
|
// convert to 3D positions and create coordMap
|
|
RDGeom::Point3DPtrVect positions3D;
|
|
for (unsigned int p = 0; p < positions.size(); ++p) {
|
|
positions3D.push_back(new RDGeom::Point3D(
|
|
(*positions[p])[0], (*positions[p])[1], (*positions[p])[2]));
|
|
}
|
|
|
|
// create the force field
|
|
ForceFields::ForceField *field;
|
|
if (useBasicKnowledge) { // ETKDG or KDG
|
|
field = DistGeom::construct3DForceField(*mmat, positions3D, bonds, angles,
|
|
expTorsionAtoms, expTorsionAngles,
|
|
improperAtoms, atomNums);
|
|
} else { // plain ETDG
|
|
field = DistGeom::constructPlain3DForceField(*mmat, positions3D, bonds,
|
|
angles, expTorsionAtoms,
|
|
expTorsionAngles, atomNums);
|
|
}
|
|
|
|
// minimize!
|
|
field->initialize();
|
|
// std::cout << "Field with torsion constraints: " << field->calcEnergy() << "
|
|
// " << ERROR_TOL << std::endl;
|
|
if (field->calcEnergy() > ERROR_TOL) {
|
|
// while (needMore) {
|
|
field->minimize(300, optimizerForceTol);
|
|
// ++nPasses;
|
|
//}
|
|
}
|
|
// std::cout << field->calcEnergy() << std::endl;
|
|
|
|
delete field;
|
|
|
|
// overwrite positions and delete the 3D ones
|
|
for (unsigned int i = 0; i < positions3D.size(); ++i) {
|
|
(*positions[i])[0] = (*positions3D[i])[0];
|
|
(*positions[i])[1] = (*positions3D[i])[1];
|
|
(*positions[i])[2] = (*positions3D[i])[2];
|
|
delete positions3D[i];
|
|
}
|
|
}
|
|
|
|
bool _embedPoints(
|
|
RDGeom::PointPtrVect *positions, const DistGeom::BoundsMatPtr mmat,
|
|
bool useRandomCoords, double boxSizeMult, bool randNegEig,
|
|
unsigned int numZeroFail, double optimizerForceTol, double basinThresh,
|
|
int seed, unsigned int maxIterations,
|
|
const DistGeom::VECT_CHIRALSET *chiralCenters, bool enforceChirality,
|
|
bool useExpTorsionAnglePrefs, bool useBasicKnowledge,
|
|
const std::vector<std::pair<int, int> > &bonds,
|
|
const std::vector<std::vector<int> > &angles,
|
|
const std::vector<std::vector<int> > &expTorsionAtoms,
|
|
const std::vector<std::pair<std::vector<int>, std::vector<double> > >
|
|
&expTorsionAngles,
|
|
const std::vector<std::vector<int> > &improperAtoms,
|
|
const std::vector<int> &atomNums) {
|
|
unsigned int nat = positions->size();
|
|
if (maxIterations == 0) {
|
|
maxIterations = 10 * nat;
|
|
}
|
|
RDNumeric::DoubleSymmMatrix distMat(nat, 0.0);
|
|
|
|
// The basin threshold just gets us into trouble when we're using
|
|
// random coordinates since it ends up ignoring 1-4 (and higher)
|
|
// interactions. This causes us to get folded-up (and self-penetrating)
|
|
// conformations for large flexible molecules
|
|
if (useRandomCoords) basinThresh = 1e8;
|
|
|
|
RDKit::double_source_type *rng = 0;
|
|
RDKit::rng_type *generator;
|
|
RDKit::uniform_double *distrib;
|
|
if (seed > 0) {
|
|
generator = new RDKit::rng_type(42u);
|
|
generator->seed(seed);
|
|
distrib = new RDKit::uniform_double(0.0, 1.0);
|
|
rng = new RDKit::double_source_type(*generator, *distrib);
|
|
} else {
|
|
rng = &RDKit::getDoubleRandomSource();
|
|
}
|
|
|
|
bool gotCoords = false;
|
|
unsigned int iter = 0;
|
|
double largestDistance = -1.0;
|
|
RDUNUSED_PARAM(largestDistance);
|
|
while ((gotCoords == false) && (iter < maxIterations)) {
|
|
++iter;
|
|
if (!useRandomCoords) {
|
|
largestDistance = DistGeom::pickRandomDistMat(*mmat, distMat, *rng);
|
|
gotCoords = DistGeom::computeInitialCoords(distMat, *positions, *rng,
|
|
randNegEig, numZeroFail);
|
|
} else {
|
|
double boxSize;
|
|
if (boxSizeMult > 0) {
|
|
boxSize = 5. * boxSizeMult;
|
|
} else {
|
|
boxSize = -1 * boxSizeMult;
|
|
}
|
|
gotCoords = DistGeom::computeRandomCoords(*positions, boxSize, *rng);
|
|
}
|
|
if (gotCoords) {
|
|
ForceFields::ForceField *field = DistGeom::constructForceField(
|
|
*mmat, *positions, *chiralCenters, 1.0, 0.1, 0, basinThresh);
|
|
unsigned int nPasses = 0;
|
|
field->initialize();
|
|
// std::cerr<<"FIELD E: "<<field->calcEnergy()<<std::endl;
|
|
if (field->calcEnergy() > ERROR_TOL) {
|
|
int needMore = 1;
|
|
while (needMore) {
|
|
needMore = field->minimize(400, optimizerForceTol);
|
|
++nPasses;
|
|
}
|
|
}
|
|
delete field;
|
|
field = NULL;
|
|
// std::cerr<<" "<<field->calcEnergy()<<" after npasses:
|
|
// "<<nPasses<<std::endl;
|
|
|
|
// Check if any of our chiral centers are badly out of whack. If so, try
|
|
// again
|
|
if (enforceChirality && chiralCenters->size() > 0) {
|
|
// check the chiral volume:
|
|
BOOST_FOREACH (DistGeom::ChiralSetPtr chiralSet, *chiralCenters) {
|
|
double vol = DistGeom::ChiralViolationContrib::calcChiralVolume(
|
|
chiralSet->d_idx1, chiralSet->d_idx2, chiralSet->d_idx3,
|
|
chiralSet->d_idx4, *positions);
|
|
double lb = chiralSet->getLowerVolumeBound();
|
|
double ub = chiralSet->getUpperVolumeBound();
|
|
if ((lb > 0 && vol < lb && (lb - vol) / lb > .2) ||
|
|
(ub < 0 && vol > ub && (vol - ub) / ub > .2)) {
|
|
// std::cerr<<" fail! ("<<chiralSet->d_idx0<<") iter: "<<iter<<"
|
|
// "<<vol<<" "<<lb<<"-"<<ub<<std::endl;
|
|
gotCoords = false;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
// now redo the minimization if we have a chiral center
|
|
// or have started from random coords. This
|
|
// time removing the chiral constraints and
|
|
// increasing the weight on the fourth dimension
|
|
if (gotCoords && (chiralCenters->size() > 0 || useRandomCoords)) {
|
|
ForceFields::ForceField *field2 = DistGeom::constructForceField(
|
|
*mmat, *positions, *chiralCenters, 0.2, 1.0, 0, basinThresh);
|
|
field2->initialize();
|
|
// std::cerr<<"FIELD2 E: "<<field2->calcEnergy()<<std::endl;
|
|
if (field2->calcEnergy() > ERROR_TOL) {
|
|
int needMore = 1;
|
|
int nPasses2 = 0;
|
|
while (needMore) {
|
|
needMore = field2->minimize(200, optimizerForceTol);
|
|
++nPasses2;
|
|
}
|
|
// std::cerr<<" "<<field2->calcEnergy()<<" after npasses2:
|
|
// "<<nPasses2<<std::endl;
|
|
}
|
|
delete field2;
|
|
}
|
|
|
|
// (ET)(K)DG
|
|
if (gotCoords && (useExpTorsionAnglePrefs || useBasicKnowledge)) {
|
|
_minimizeWithExpTorsions(*positions, mmat, optimizerForceTol,
|
|
basinThresh, bonds, angles, expTorsionAtoms,
|
|
expTorsionAngles, improperAtoms, atomNums,
|
|
useBasicKnowledge);
|
|
}
|
|
|
|
// test if chirality is correct
|
|
if (enforceChirality && gotCoords && (chiralCenters->size() > 0)) {
|
|
// "distance matrix" chirality test
|
|
std::set<int> atoms;
|
|
BOOST_FOREACH (DistGeom::ChiralSetPtr chiralSet, *chiralCenters) {
|
|
if (chiralSet->d_idx0 != chiralSet->d_idx4) {
|
|
atoms.insert(chiralSet->d_idx0);
|
|
atoms.insert(chiralSet->d_idx1);
|
|
atoms.insert(chiralSet->d_idx2);
|
|
atoms.insert(chiralSet->d_idx3);
|
|
atoms.insert(chiralSet->d_idx4);
|
|
}
|
|
}
|
|
std::vector<int> atomsToCheck(atoms.begin(), atoms.end());
|
|
if (atomsToCheck.size() > 0) {
|
|
if (!_boundsFulfilled(atomsToCheck, *mmat, *positions)) {
|
|
gotCoords = false;
|
|
}
|
|
}
|
|
|
|
// "center in volume" chirality test
|
|
if (gotCoords) {
|
|
BOOST_FOREACH (DistGeom::ChiralSetPtr chiralSet, *chiralCenters) {
|
|
// it could happen that the centroid is outside the volume defined
|
|
// by the other
|
|
// four points. That is also a fail.
|
|
if (!_centerInVolume(chiralSet, *positions)) {
|
|
// std::cerr<<" fail2! ("<<chiralSet->d_idx0<<") iter:
|
|
// "<<iter<<std::endl;
|
|
gotCoords = false;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
} // if(gotCoords)
|
|
} // while
|
|
if (seed > 0 && rng) {
|
|
delete rng;
|
|
delete generator;
|
|
delete distrib;
|
|
}
|
|
return gotCoords;
|
|
}
|
|
|
|
void _findChiralSets(const ROMol &mol,
|
|
DistGeom::VECT_CHIRALSET &chiralCenters) {
|
|
ROMol::ConstAtomIterator ati;
|
|
INT_VECT nbrs;
|
|
ROMol::OEDGE_ITER beg, end;
|
|
// Atom *oatom;
|
|
for (ati = mol.beginAtoms(); ati != mol.endAtoms(); ati++) {
|
|
if ((*ati)->getAtomicNum() != 1) { // skip hydrogens
|
|
Atom::ChiralType chiralType = (*ati)->getChiralTag();
|
|
if (chiralType == Atom::CHI_TETRAHEDRAL_CW ||
|
|
chiralType == Atom::CHI_TETRAHEDRAL_CCW) {
|
|
// make a chiral set from the neighbors
|
|
nbrs.clear();
|
|
nbrs.reserve(4);
|
|
// find the neighbors of this atom and enter them into the
|
|
// nbr list
|
|
boost::tie(beg, end) = mol.getAtomBonds(*ati);
|
|
while (beg != end) {
|
|
nbrs.push_back(mol[*beg]->getOtherAtom(*ati)->getIdx());
|
|
++beg;
|
|
}
|
|
// if we have less than 4 heavy atoms as neighbors,
|
|
// we need to include the chiral center into the mix
|
|
// we should at least have 3 though
|
|
bool includeSelf = false;
|
|
RDUNUSED_PARAM(includeSelf);
|
|
CHECK_INVARIANT(nbrs.size() >= 3, "Cannot be a chiral center");
|
|
|
|
if (nbrs.size() < 4) {
|
|
nbrs.insert(nbrs.end(), (*ati)->getIdx());
|
|
includeSelf = true;
|
|
}
|
|
|
|
// now create a chiral set and set the upper and lower bound on the
|
|
// volume
|
|
if (chiralType == Atom::CHI_TETRAHEDRAL_CCW) {
|
|
// postive chiral volume
|
|
DistGeom::ChiralSet *cset = new DistGeom::ChiralSet(
|
|
(*ati)->getIdx(), nbrs[0], nbrs[1], nbrs[2], nbrs[3], 5.0, 100.0);
|
|
DistGeom::ChiralSetPtr cptr(cset);
|
|
chiralCenters.push_back(cptr);
|
|
} else {
|
|
DistGeom::ChiralSet *cset =
|
|
new DistGeom::ChiralSet((*ati)->getIdx(), nbrs[0], nbrs[1],
|
|
nbrs[2], nbrs[3], -100.0, -5.0);
|
|
DistGeom::ChiralSetPtr cptr(cset);
|
|
chiralCenters.push_back(cptr);
|
|
}
|
|
} // if block -chirality check
|
|
} // if block - heavy atom check
|
|
} // for loop over atoms
|
|
} // end of _findChiralSets
|
|
|
|
void _fillAtomPositions(RDGeom::Point3DConstPtrVect &pts,
|
|
const Conformer &conf) {
|
|
unsigned int na = conf.getNumAtoms();
|
|
pts.clear();
|
|
unsigned int ai;
|
|
pts.reserve(na);
|
|
for (ai = 0; ai < na; ++ai) {
|
|
pts.push_back(&conf.getAtomPos(ai));
|
|
}
|
|
}
|
|
|
|
bool _isConfFarFromRest(const ROMol &mol, const Conformer &conf,
|
|
double threshold) {
|
|
// NOTE: it is tempting to use some triangle inequality to prune
|
|
// conformations here but some basic testing has shown very
|
|
// little advantage and given that the time for pruning fades in
|
|
// comparison to embedding - we will use a simple for loop below
|
|
// over all conformation until we find a match
|
|
ROMol::ConstConformerIterator confi;
|
|
|
|
RDGeom::Point3DConstPtrVect refPoints, prbPoints;
|
|
_fillAtomPositions(refPoints, conf);
|
|
|
|
bool res = true;
|
|
unsigned int na = conf.getNumAtoms();
|
|
double ssrThres = na * threshold * threshold;
|
|
|
|
RDGeom::Transform3D trans;
|
|
double ssr;
|
|
for (confi = mol.beginConformers(); confi != mol.endConformers(); confi++) {
|
|
_fillAtomPositions(prbPoints, *(*confi));
|
|
ssr = RDNumeric::Alignments::AlignPoints(refPoints, prbPoints, trans);
|
|
if (ssr < ssrThres) {
|
|
res = false;
|
|
break;
|
|
}
|
|
}
|
|
return res;
|
|
}
|
|
|
|
int EmbedMolecule(ROMol &mol, unsigned int maxIterations, int seed,
|
|
bool clearConfs, bool useRandomCoords, double boxSizeMult,
|
|
bool randNegEig, unsigned int numZeroFail,
|
|
const std::map<int, RDGeom::Point3D> *coordMap,
|
|
double optimizerForceTol, bool ignoreSmoothingFailures,
|
|
bool enforceChirality, bool useExpTorsionAnglePrefs,
|
|
bool useBasicKnowledge, bool verbose, double basinThresh) {
|
|
INT_VECT confIds;
|
|
EmbedMultipleConfs(mol, confIds, 1, 1, maxIterations, seed, clearConfs,
|
|
useRandomCoords, boxSizeMult, randNegEig, numZeroFail,
|
|
-1.0, coordMap, optimizerForceTol, ignoreSmoothingFailures,
|
|
enforceChirality, useExpTorsionAnglePrefs,
|
|
useBasicKnowledge, verbose, basinThresh);
|
|
|
|
int res;
|
|
if (confIds.size()) {
|
|
res = confIds[0];
|
|
} else {
|
|
res = -1;
|
|
}
|
|
return res;
|
|
}
|
|
|
|
void adjustBoundsMatFromCoordMap(
|
|
DistGeom::BoundsMatPtr mmat, unsigned int nAtoms,
|
|
const std::map<int, RDGeom::Point3D> *coordMap) {
|
|
RDUNUSED_PARAM(nAtoms);
|
|
// std::cerr<<std::endl;
|
|
// for(unsigned int i=0;i<nAtoms;++i){
|
|
// for(unsigned int j=0;j<nAtoms;++j){
|
|
// std::cerr<<" "<<std::setprecision(3)<<mmat->getVal(i,j);
|
|
// }
|
|
// std::cerr<<std::endl;
|
|
// }
|
|
// std::cerr<<std::endl;
|
|
for (std::map<int, RDGeom::Point3D>::const_iterator iIt = coordMap->begin();
|
|
iIt != coordMap->end(); ++iIt) {
|
|
int iIdx = iIt->first;
|
|
const RDGeom::Point3D &iPoint = iIt->second;
|
|
|
|
std::map<int, RDGeom::Point3D>::const_iterator jIt = iIt;
|
|
while (++jIt != coordMap->end()) {
|
|
int jIdx = jIt->first;
|
|
const RDGeom::Point3D &jPoint = jIt->second;
|
|
double dist = (iPoint - jPoint).length();
|
|
mmat->setUpperBound(iIdx, jIdx, dist);
|
|
mmat->setLowerBound(iIdx, jIdx, dist);
|
|
}
|
|
}
|
|
// std::cerr<<std::endl;
|
|
// for(unsigned int i=0;i<nAtoms;++i){
|
|
// for(unsigned int j=0;j<nAtoms;++j){
|
|
// std::cerr<<" "<<std::setprecision(3)<<mmat->getVal(i,j);
|
|
// }
|
|
// std::cerr<<std::endl;
|
|
// }
|
|
// std::cerr<<std::endl;
|
|
}
|
|
|
|
namespace detail {
|
|
typedef struct {
|
|
boost::dynamic_bitset<> *confsOk;
|
|
bool fourD;
|
|
INT_VECT *fragMapping;
|
|
std::vector<Conformer *> *confs;
|
|
unsigned int fragIdx;
|
|
DistGeom::BoundsMatPtr mmat;
|
|
bool useRandomCoords;
|
|
double boxSizeMult;
|
|
bool randNegEig;
|
|
unsigned int numZeroFail;
|
|
double optimizerForceTol;
|
|
double basinThresh;
|
|
int seed;
|
|
unsigned int maxIterations;
|
|
DistGeom::VECT_CHIRALSET const *chiralCenters;
|
|
bool enforceChirality;
|
|
bool useExpTorsionAnglePrefs;
|
|
bool useBasicKnowledge;
|
|
std::vector<std::pair<int, int> > *bonds;
|
|
std::vector<std::vector<int> > *angles;
|
|
std::vector<std::vector<int> > *expTorsionAtoms;
|
|
std::vector<std::pair<std::vector<int>, std::vector<double> > >
|
|
*expTorsionAngles;
|
|
std::vector<std::vector<int> > *improperAtoms;
|
|
std::vector<int> *atomNums;
|
|
} EmbedArgs;
|
|
void embedHelper_(int threadId, int numThreads, EmbedArgs *eargs) {
|
|
unsigned int nAtoms = eargs->mmat->numRows();
|
|
RDGeom::PointPtrVect positions;
|
|
for (unsigned int i = 0; i < nAtoms; ++i) {
|
|
if (eargs->fourD) {
|
|
positions.push_back(new RDGeom::PointND(4));
|
|
} else {
|
|
positions.push_back(new RDGeom::Point3D());
|
|
}
|
|
}
|
|
for (size_t ci = 0; ci < eargs->confs->size(); ci++) {
|
|
if (rdcast<int>(ci % numThreads) != threadId) continue;
|
|
if (!(*eargs->confsOk)[ci]) {
|
|
// if one of the fragments here has already failed, there's no
|
|
// sense in embedding this one
|
|
continue;
|
|
}
|
|
bool gotCoords = _embedPoints(
|
|
&positions, eargs->mmat, eargs->useRandomCoords, eargs->boxSizeMult,
|
|
eargs->randNegEig, eargs->numZeroFail, eargs->optimizerForceTol,
|
|
eargs->basinThresh, (ci + 1) * eargs->seed, eargs->maxIterations,
|
|
eargs->chiralCenters, eargs->enforceChirality,
|
|
eargs->useExpTorsionAnglePrefs, eargs->useBasicKnowledge, *eargs->bonds,
|
|
*eargs->angles, *eargs->expTorsionAtoms, *eargs->expTorsionAngles,
|
|
*eargs->improperAtoms, *eargs->atomNums);
|
|
|
|
if (gotCoords) {
|
|
Conformer *conf = (*eargs->confs)[ci];
|
|
unsigned int fragAtomIdx = 0;
|
|
for (unsigned int i = 0; i < (*eargs->confs)[0]->getNumAtoms(); ++i) {
|
|
if ((*eargs->fragMapping)[i] == static_cast<int>(eargs->fragIdx)) {
|
|
conf->setAtomPos(i, RDGeom::Point3D((*positions[fragAtomIdx])[0],
|
|
(*positions[fragAtomIdx])[1],
|
|
(*positions[fragAtomIdx])[2]));
|
|
++fragAtomIdx;
|
|
}
|
|
}
|
|
} else {
|
|
(*eargs->confsOk)[ci] = 0;
|
|
}
|
|
}
|
|
for (unsigned int i = 0; i < nAtoms; ++i) {
|
|
delete positions[i];
|
|
}
|
|
}
|
|
} // end of namespace detail
|
|
|
|
void EmbedMultipleConfs(ROMol &mol, INT_VECT &res, unsigned int numConfs,
|
|
int numThreads, unsigned int maxIterations, int seed,
|
|
bool clearConfs, bool useRandomCoords,
|
|
double boxSizeMult, bool randNegEig,
|
|
unsigned int numZeroFail, double pruneRmsThresh,
|
|
const std::map<int, RDGeom::Point3D> *coordMap,
|
|
double optimizerForceTol, bool ignoreSmoothingFailures,
|
|
bool enforceChirality, bool useExpTorsionAnglePrefs,
|
|
bool useBasicKnowledge, bool verbose,
|
|
double basinThresh) {
|
|
if (!mol.getNumAtoms()) {
|
|
throw ValueErrorException("molecule has no atoms");
|
|
}
|
|
|
|
INT_VECT fragMapping;
|
|
std::vector<ROMOL_SPTR> molFrags =
|
|
MolOps::getMolFrags(mol, true, &fragMapping);
|
|
if (molFrags.size() > 1 && coordMap) {
|
|
BOOST_LOG(rdWarningLog)
|
|
<< "Constrained conformer generation (via the coordMap argument) does "
|
|
"not work with molecules that have multiple fragments."
|
|
<< std::endl;
|
|
coordMap = 0;
|
|
}
|
|
std::vector<Conformer *> confs;
|
|
confs.reserve(numConfs);
|
|
for (unsigned int i = 0; i < numConfs; ++i) {
|
|
confs.push_back(new Conformer(mol.getNumAtoms()));
|
|
}
|
|
boost::dynamic_bitset<> confsOk(numConfs);
|
|
confsOk.set();
|
|
|
|
if (clearConfs) {
|
|
res.clear();
|
|
mol.clearConformers();
|
|
}
|
|
|
|
for (unsigned int fragIdx = 0; fragIdx < molFrags.size(); ++fragIdx) {
|
|
ROMOL_SPTR piece = molFrags[fragIdx];
|
|
unsigned int nAtoms = piece->getNumAtoms();
|
|
DistGeom::BoundsMatrix *mat = new DistGeom::BoundsMatrix(nAtoms);
|
|
DistGeom::BoundsMatPtr mmat(mat);
|
|
initBoundsMat(mmat);
|
|
|
|
double tol = 0.0;
|
|
std::vector<std::vector<int> > expTorsionAtoms;
|
|
std::vector<std::pair<std::vector<int>, std::vector<double> > >
|
|
expTorsionAngles;
|
|
std::vector<std::vector<int> > improperAtoms;
|
|
std::vector<std::pair<int, int> > bonds;
|
|
std::vector<std::vector<int> > angles;
|
|
std::vector<int> atomNums(nAtoms);
|
|
if (useExpTorsionAnglePrefs || useBasicKnowledge) {
|
|
ForceFields::CrystalFF::getExperimentalTorsions(
|
|
*piece, expTorsionAtoms, expTorsionAngles, improperAtoms,
|
|
useExpTorsionAnglePrefs, useBasicKnowledge, verbose);
|
|
setTopolBounds(*piece, mmat, bonds, angles, true, false);
|
|
for (unsigned int i = 0; i < nAtoms; ++i) {
|
|
atomNums[i] = (*piece).getAtomWithIdx(i)->getAtomicNum();
|
|
}
|
|
} else {
|
|
setTopolBounds(*piece, mmat, true, false);
|
|
}
|
|
if (coordMap) {
|
|
adjustBoundsMatFromCoordMap(mmat, nAtoms, coordMap);
|
|
tol = 0.05;
|
|
}
|
|
if (!DistGeom::triangleSmoothBounds(mmat, tol)) {
|
|
// ok this bound matrix failed to triangle smooth - re-compute the bounds
|
|
// matrix
|
|
// without 15 bounds and with VDW scaling
|
|
initBoundsMat(mmat);
|
|
setTopolBounds(*piece, mmat, false, true);
|
|
|
|
if (coordMap) {
|
|
adjustBoundsMatFromCoordMap(mmat, nAtoms, coordMap);
|
|
}
|
|
|
|
// try triangle smoothing again
|
|
if (!DistGeom::triangleSmoothBounds(mmat, tol)) {
|
|
// ok, we're not going to be able to smooth this,
|
|
if (ignoreSmoothingFailures) {
|
|
// proceed anyway with the more relaxed bounds matrix
|
|
initBoundsMat(mmat);
|
|
setTopolBounds(*piece, mmat, false, true);
|
|
|
|
if (coordMap) {
|
|
adjustBoundsMatFromCoordMap(mmat, nAtoms, coordMap);
|
|
}
|
|
} else {
|
|
BOOST_LOG(rdWarningLog)
|
|
<< "Could not triangle bounds smooth molecule." << std::endl;
|
|
return;
|
|
}
|
|
}
|
|
}
|
|
#if 0
|
|
for(unsigned int li=0;li<piece->getNumAtoms();++li){
|
|
for(unsigned int lj=li+1;lj<piece->getNumAtoms();++lj){
|
|
std::cerr<<" ("<<li<<","<<lj<<"): "<<mat->getLowerBound(li,lj)<<" -> "<<mat->getUpperBound(li,lj)<<std::endl;
|
|
}
|
|
}
|
|
#endif
|
|
// find all the chiral centers in the molecule
|
|
DistGeom::VECT_CHIRALSET chiralCenters;
|
|
MolOps::assignStereochemistry(*piece);
|
|
_findChiralSets(*piece, chiralCenters);
|
|
|
|
// if we have any chiral centers or are using random coordinates, we will
|
|
// first embed the molecule in four dimensions, otherwise we will use 3D
|
|
bool fourD = false;
|
|
if (useRandomCoords || chiralCenters.size() > 0) {
|
|
fourD = true;
|
|
}
|
|
#ifdef RDK_THREADSAFE_SSS
|
|
boost::thread_group tg;
|
|
#endif
|
|
numThreads = getNumThreadsToUse(numThreads);
|
|
|
|
detail::EmbedArgs eargs = {&confsOk,
|
|
fourD,
|
|
&fragMapping,
|
|
&confs,
|
|
fragIdx,
|
|
mmat,
|
|
useRandomCoords,
|
|
boxSizeMult,
|
|
randNegEig,
|
|
numZeroFail,
|
|
optimizerForceTol,
|
|
basinThresh,
|
|
seed,
|
|
maxIterations,
|
|
&chiralCenters,
|
|
enforceChirality,
|
|
useExpTorsionAnglePrefs,
|
|
useBasicKnowledge,
|
|
&bonds,
|
|
&angles,
|
|
&expTorsionAtoms,
|
|
&expTorsionAngles,
|
|
&improperAtoms,
|
|
&atomNums};
|
|
if (numThreads == 1) {
|
|
detail::embedHelper_(0, 1, &eargs);
|
|
}
|
|
#ifdef RDK_THREADSAFE_SSS
|
|
else {
|
|
for (int tid = 0; tid < numThreads; ++tid) {
|
|
tg.add_thread(
|
|
new boost::thread(detail::embedHelper_, tid, numThreads, &eargs));
|
|
}
|
|
tg.join_all();
|
|
}
|
|
#endif
|
|
}
|
|
for (unsigned int ci = 0; ci < confs.size(); ++ci) {
|
|
Conformer *conf = confs[ci];
|
|
if (confsOk[ci]) {
|
|
// check if we are pruning away conformations and
|
|
// a closeby conformation has already been chosen :
|
|
if (pruneRmsThresh > 0.0 &&
|
|
!_isConfFarFromRest(mol, *conf, pruneRmsThresh)) {
|
|
delete conf;
|
|
} else {
|
|
int confId = (int)mol.addConformer(conf, true);
|
|
res.push_back(confId);
|
|
}
|
|
} else {
|
|
delete conf;
|
|
}
|
|
}
|
|
}
|
|
|
|
INT_VECT EmbedMultipleConfs(
|
|
ROMol &mol, unsigned int numConfs, unsigned int maxIterations, int seed,
|
|
bool clearConfs, bool useRandomCoords, double boxSizeMult, bool randNegEig,
|
|
unsigned int numZeroFail, double pruneRmsThresh,
|
|
const std::map<int, RDGeom::Point3D> *coordMap, double optimizerForceTol,
|
|
bool ignoreSmoothingFailures, bool enforceChirality,
|
|
bool useExpTorsionAnglePrefs, bool useBasicKnowledge, bool verbose,
|
|
double basinThresh) {
|
|
INT_VECT res;
|
|
EmbedMultipleConfs(
|
|
mol, res, numConfs, 1, maxIterations, seed, clearConfs, useRandomCoords,
|
|
boxSizeMult, randNegEig, numZeroFail, pruneRmsThresh, coordMap,
|
|
optimizerForceTol, ignoreSmoothingFailures, enforceChirality,
|
|
useExpTorsionAnglePrefs, useBasicKnowledge, verbose, basinThresh);
|
|
return res;
|
|
}
|
|
} // end of namespace DGeomHelpers
|
|
} // end of namespace RDKit
|