pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-03 21:44:30 +08:00
Code Issues Packages Projects Releases Wiki Activity

Refactor distgeom minimizations (#7652)

Browse Source 
* refactor skeleton

* fix contribs

* adapt documentation

* clarify function names

* make variable names more readable and update documentation

* switch boost dynamic bitset with std::vector<bool>

* clearer const name

* fix typo and go back to dynamic bitset

* go back to old function signature

* add inline overload with fewer arguments

* remove unecessary doc lines and add tol to distance directly in function call

* use constexpr for known distance force constant

* update documentation

* enhance documentation

* english

Co-authored-by: Greg Landrum <greg.landrum@gmail.com>

* make changes requested by @greglandrum

* make all constexprs

---------

Co-authored-by: Greg Landrum <greg.landrum@gmail.com>
This commit is contained in:
nmaeder
2024-07-25 06:14:12 +02:00
committed by GitHub
parent bd49c8ae81
commit 3f3ff3858a
3 changed files with 315 additions and 289 deletions
Download Patch File Download Diff File Expand all files Collapse all files

570
Code/DistGeom/DistGeomUtils.cpp
View File

@@ -28,7 +28,9 @@
#include <ForceField/MMFF/Nonbonded.h>
namespace DistGeom {
const double EIGVAL_TOL = 0.001;
constexpr double EIGVAL_TOL = 0.001;
constexpr double KNOWN_DIST_TOL = 0.01;
constexpr double KNOWN_DIST_FORCE_CONSTANT = 100;
double pickRandomDistMat(const BoundsMatrix &mmat,
RDNumeric::SymmMatrix<double> &distMat, int seed) {
@@ -187,9 +189,8 @@ ForceFields::ForceField *constructForceField(
unsigned int N = mmat.numRows();
CHECK_INVARIANT(N == positions.size(), "");
auto *field = new ForceFields::ForceField(positions[0]->dimension());
for (unsigned int i = 0; i < N; i++) {
field->positions().push_back(positions[i]);
}
field->positions().insert(field->positions().begin(), positions.begin(),
positions.end());
auto contrib = new DistViolationContribs(field);
for (unsigned int i = 1; i < N; i++) {
@@ -251,278 +252,23 @@ ForceFields::ForceField *constructForceField(
return field;
} // constructForceField
ForceFields::ForceField *construct3DForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails) {
unsigned int N = mmat.numRows();
CHECK_INVARIANT(N == positions.size(), "");
CHECK_INVARIANT(etkdgDetails.expTorsionAtoms.size() ==
etkdgDetails.expTorsionAngles.size(),
"");
auto *field = new ForceFields::ForceField(positions[0]->dimension());
field->positions().insert(field->positions().begin(), positions.begin(),
positions.end());
//! Add basic knowledge improper torsion contributions to a force field
/*!
// keep track which atoms are 1,2- or 1,3-restrained
boost::dynamic_bitset<> atomPairs(N * N);
\param ff Force field to add contributions to
\param forceScalingFactor Force scaling factor to use in inversion contrib
\param improperAtoms Indices of atoms to be used in improper torsion
terms.
\param isImproperConstrained bit vector with length of total num atoms of
the molecule where index of every central atom of improper torsion is set to
one
// don't add 1-3 Distances constraints for angles where the
// central atom of the angle is the central atom of an improper torsion.
boost::dynamic_bitset<> dont13Constrain(N);
// torsion constraints
for (unsigned int t = 0; t < etkdgDetails.expTorsionAtoms.size(); ++t) {
int i = etkdgDetails.expTorsionAtoms[t][0];
int j = etkdgDetails.expTorsionAtoms[t][1];
int k = etkdgDetails.expTorsionAtoms[t][2];
int l = etkdgDetails.expTorsionAtoms[t][3];
if (i < l) {
atomPairs[i * N + l] = 1;
} else {
atomPairs[l * N + i] = 1;
}
// etkdgDetails.expTorsionAngles[t][0] = (signs, V's)
auto *contrib = new ForceFields::CrystalFF::TorsionAngleContribM6(
field, i, j, k, l, etkdgDetails.expTorsionAngles[t].second,
etkdgDetails.expTorsionAngles[t].first);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
} // torsion constraints
// improper torsions / out-of-plane bend / inversion
double oobForceScalingFactor = 10.0;
for (const auto &improperAtom : etkdgDetails.improperAtoms) {
std::vector<int> n(4);
for (unsigned int i = 0; i < 3; ++i) {
n[1] = 1;
switch (i) {
case 0:
n[0] = 0;
n[2] = 2;
n[3] = 3;
break;
case 1:
n[0] = 0;
n[2] = 3;
n[3] = 2;
break;
case 2:
n[0] = 2;
n[2] = 3;
n[3] = 0;
break;
}
auto *contrib = new ForceFields::UFF::InversionContrib(
field, improperAtom[n[0]], improperAtom[n[1]], improperAtom[n[2]],
improperAtom[n[3]], improperAtom[4],
static_cast<bool>(improperAtom[5]), oobForceScalingFactor);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
dont13Constrain[improperAtom[n[1]]] = 1;
}
}
constexpr double knownDistanceConstraintForce = 100.0;
double fdist = knownDistanceConstraintForce; // force constant
// 1,2 distance constraints
for (const auto &bnd : etkdgDetails.bonds) {
unsigned int i = bnd.first;
unsigned int j = bnd.second;
if (i < j) {
atomPairs[i * N + j] = 1;
} else {
atomPairs[j * N + i] = 1;
}
double d = ((*positions[i]) - (*positions[j])).length();
double l = d - 0.01;
double u = d + 0.01;
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
field, i, j, l, u, fdist);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
// 1,3 distance constraints
for (const auto &angle : etkdgDetails.angles) {
unsigned int i = angle[0];
unsigned int j = angle[1];
unsigned int k = angle[2];
if (i < k) {
atomPairs[i * N + k] = 1;
} else {
atomPairs[k * N + i] = 1;
}
// check for triple bonds
if (angle[3]) {
auto *contrib = new ForceFields::UFF::AngleConstraintContrib(
field, i, j, k, 179.0, 180.0, 1);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
} else if (!dont13Constrain.test(j)) {
double d = ((*positions[i]) - (*positions[k])).length();
double l = d - 0.01;
double u = d + 0.01;
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
field, i, k, l, u, fdist);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
}
// minimum distance for all other atom pairs that aren't constrained
for (unsigned int i = 1; i < N; ++i) {
for (unsigned int j = 0; j < i; ++j) {
if (!atomPairs[j * N + i]) {
fdist = etkdgDetails.boundsMatForceScaling * 10.0;
double l = mmat.getLowerBound(i, j);
double u = mmat.getUpperBound(i, j);
if (!etkdgDetails.constrainedAtoms.empty() &&
etkdgDetails.constrainedAtoms[i] &&
etkdgDetails.constrainedAtoms[j]) {
// we're constrained, so use very tight bounds
l = u = ((*positions[i]) - (*positions[j])).length();
constexpr double INCR = 0.01;
l -= INCR;
u += INCR;
fdist = knownDistanceConstraintForce;
}
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
field, i, j, l, u, fdist);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
}
}
return field;
} // construct3DForceField
ForceFields::ForceField *construct3DForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails,
const std::map<std::pair<unsigned int, unsigned int>, double> &CPCI) {
auto *field = construct3DForceField(mmat, positions, etkdgDetails);
bool is1_4 = false;
// double dielConst = 1.0;
boost::uint8_t dielModel = 1;
for (const auto &charge : CPCI) {
auto *contrib = new ForceFields::MMFF::EleContrib(
field, charge.first.first, charge.first.second, charge.second,
dielModel, is1_4);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
return field;
}
ForceFields::ForceField *constructPlain3DForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails) {
unsigned int N = mmat.numRows();
CHECK_INVARIANT(N == positions.size(), "");
CHECK_INVARIANT(etkdgDetails.expTorsionAtoms.size() ==
etkdgDetails.expTorsionAngles.size(),
"");
auto *field = new ForceFields::ForceField(positions[0]->dimension());
field->positions().insert(field->positions().begin(), positions.begin(),
positions.end());
// keep track which atoms are 1,2- or 1,3-restrained
boost::dynamic_bitset<> atomPairs(N * N);
// torsion constraints
for (unsigned int t = 0; t < etkdgDetails.expTorsionAtoms.size(); ++t) {
int i = etkdgDetails.expTorsionAtoms[t][0];
int j = etkdgDetails.expTorsionAtoms[t][1];
int k = etkdgDetails.expTorsionAtoms[t][2];
int l = etkdgDetails.expTorsionAtoms[t][3];
if (i < l) {
atomPairs[i * N + l] = 1;
} else {
atomPairs[l * N + i] = 1;
}
// etkdgDetails.expTorsionAngles[t][0] = (signs, V's)
auto *contrib = new ForceFields::CrystalFF::TorsionAngleContribM6(
field, i, j, k, l, etkdgDetails.expTorsionAngles[t].second,
etkdgDetails.expTorsionAngles[t].first);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
} // torsion constraints
constexpr double knownDistanceConstraintForce = 100.0;
double fdist = knownDistanceConstraintForce; // force constant
// 1,2 distance constraints
for (const auto &bnd : etkdgDetails.bonds) {
unsigned int i = bnd.first;
unsigned int j = bnd.second;
if (i < j) {
atomPairs[i * N + j] = 1;
} else {
atomPairs[j * N + i] = 1;
}
double d = ((*positions[i]) - (*positions[j])).length();
double l = d - 0.01;
double u = d + 0.01;
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
field, i, j, l, u, fdist);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
// 1,3 distance constraints
for (const auto &angle : etkdgDetails.angles) {
unsigned int i = angle[0];
unsigned int k = angle[2];
if (i < k) {
atomPairs[i * N + k] = 1;
} else {
atomPairs[k * N + i] = 1;
}
double d = ((*positions[i]) - (*positions[k])).length();
double l = d - 0.01;
double u = d + 0.01;
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
field, i, k, l, u, fdist);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
// minimum distance for all other atom pairs
for (unsigned int i = 1; i < N; ++i) {
for (unsigned int j = 0; j < i; ++j) {
if (!atomPairs[j * N + i]) {
fdist = etkdgDetails.boundsMatForceScaling * 10.0;
double l = mmat.getLowerBound(i, j);
double u = mmat.getUpperBound(i, j);
if (!etkdgDetails.constrainedAtoms.empty() &&
etkdgDetails.constrainedAtoms[i] &&
etkdgDetails.constrainedAtoms[j]) {
// we're constrained, so use very tight bounds
l = u = ((*positions[i]) - (*positions[j])).length();
constexpr double INCR = 0.01;
l -= INCR;
u += INCR;
fdist = knownDistanceConstraintForce;
}
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
field, i, j, l, u, fdist);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
}
}
return field;
} // constructPlain3DForceField
ForceFields::ForceField *construct3DImproperForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const std::vector<std::vector<int>> &improperAtoms,
const std::vector<std::vector<int>> &angles,
const std::vector<int> &atomNums) {
RDUNUSED_PARAM(atomNums);
unsigned int N = mmat.numRows();
CHECK_INVARIANT(N == positions.size(), "");
auto *field = new ForceFields::ForceField(positions[0]->dimension());
field->positions().insert(field->positions().begin(), positions.begin(),
positions.end());
// improper torsions / out-of-plane bend / inversion
double oobForceScalingFactor = 10.0;
*/
void addImproperTorsionTerms(ForceFields::ForceField *ff,
double forceScalingFactor,
const std::vector<std::vector<int>> &improperAtoms,
boost::dynamic_bitset<> &isImproperConstrained) {
PRECONDITION(ff, "bad force field");
for (const auto &improperAtom : improperAtoms) {
std::vector<int> n(4);
for (unsigned int i = 0; i < 3; ++i) {
@@ -546,14 +292,277 @@ ForceFields::ForceField *construct3DImproperForceField(
n[3] = 0;
break;
}
auto *contrib = new ForceFields::UFF::InversionContrib(
field, improperAtom[n[0]], improperAtom[n[1]], improperAtom[n[2]],
improperAtom[n[3]], improperAtom[4],
static_cast<bool>(improperAtom[5]), oobForceScalingFactor);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
auto *contrib = new ForceFields::UFF::InversionContrib(
ff, improperAtom[n[0]], improperAtom[n[1]], improperAtom[n[2]],
improperAtom[n[3]], improperAtom[4],
static_cast<bool>(improperAtom[5]), forceScalingFactor);
ff->contribs().emplace_back(contrib);
isImproperConstrained[improperAtom[n[1]]] = 1;
}
}
}
//! Add experimental torsion angle contributions to a force field
/*!
\param ff Force field to add contributions to
\param etkdgDetails Contains information about the ETKDG force field
\param atomPairs bit set for every atom pair in the molecule where
a bit is set to one when the atom pair are the end atoms of a torsion
angle contribution
\param numAtoms number of atoms in the molecule
*/
void addExperimentalTorsionTerms(
ForceFields::ForceField *ff,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails,
boost::dynamic_bitset<> &atomPairs, unsigned int numAtoms) {
PRECONDITION(ff, "bad force field");
for (unsigned int t = 0; t < etkdgDetails.expTorsionAtoms.size(); ++t) {
int i = etkdgDetails.expTorsionAtoms[t][0];
int j = etkdgDetails.expTorsionAtoms[t][1];
int k = etkdgDetails.expTorsionAtoms[t][2];
int l = etkdgDetails.expTorsionAtoms[t][3];
if (i < l) {
atomPairs[i * numAtoms + l] = 1;
} else {
atomPairs[l * numAtoms + i] = 1;
}
auto *contrib = new ForceFields::CrystalFF::TorsionAngleContribM6(
ff, i, j, k, l, etkdgDetails.expTorsionAngles[t].second,
etkdgDetails.expTorsionAngles[t].first);
ff->contribs().emplace_back(contrib);
}
}
//! Add bond constraints with padding at current positions to force field
/*!
\param ff Force field to add contributions to
\param etkdgDetails Contains information about the ETKDG force field
\param atomPairs bit set for every atom pair in the molecule where
a bit is set to one when the atom pair is a bond that is constrained here
\param positions A vector of pointers to 3D Points to write out the
resulting coordinates
\param forceConstant force constant with which to constrain bond distances
\param numAtoms number of atoms in molecule
*/
void add12Terms(ForceFields::ForceField *ff,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails,
boost::dynamic_bitset<> &atomPairs,
RDGeom::Point3DPtrVect &positions, double forceConstant,
unsigned int numAtoms) {
PRECONDITION(ff, "bad force field");
for (const auto &bond : etkdgDetails.bonds) {
unsigned int i = bond.first;
unsigned int j = bond.second;
if (i < j) {
atomPairs[i * numAtoms + j] = 1;
} else {
atomPairs[j * numAtoms + i] = 1;
}
double d = ((*positions[i]) - (*positions[j])).length();
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
ff, i, j, d - KNOWN_DIST_TOL, d + KNOWN_DIST_TOL, forceConstant);
ff->contribs().emplace_back(contrib);
}
}
//! Add 1-3 distance constraints with padding at current positions to force
/// field
/*!
\param ff Force field to add contributions to
\param etkdgDetails Contains information about the ETKDG force field
\param atomPairs bit set for every atom pair in the molecule where
a bit is set to one when the atom pair is the both end atoms of a 13
contribution that is constrained here
\param positions A vector of pointers to 3D Points to write out the resulting
coordinates \param forceConstant force constant with which to constrain bond
distances \param isImproperConstrained bit vector with length of total num
atoms of the molecule where index of every central atom of improper torsion is
set to one \param useBasicKnowledge whether to use basic knowledge terms
\param numAtoms number of atoms in molecule
*/
void add13Terms(ForceFields::ForceField *ff,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails,
boost::dynamic_bitset<> &atomPairs,
RDGeom::Point3DPtrVect &positions, double forceConstant,
const boost::dynamic_bitset<> &isImproperConstrained,
bool useBasicKnowledge, unsigned int numAtoms) {
PRECONDITION(ff, "bad force field");
for (const auto &angle : etkdgDetails.angles) {
unsigned int i = angle[0];
unsigned int j = angle[1];
unsigned int k = angle[2];
if (i < k) {
atomPairs[i * numAtoms + k] = 1;
} else {
atomPairs[k * numAtoms + i] = 1;
}
// check for triple bonds
if (useBasicKnowledge && angle[3]) {
auto *contrib = new ForceFields::UFF::AngleConstraintContrib(
ff, i, j, k, 179.0, 180.0, 1);
ff->contribs().emplace_back(contrib);
} else if (!isImproperConstrained[j]) {
double d = ((*positions[i]) - (*positions[k])).length();
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
ff, i, k, d - KNOWN_DIST_TOL, d + KNOWN_DIST_TOL, forceConstant);
ff->contribs().emplace_back(contrib);
}
}
}
//! Add long distance constraints to bounds matrix borders or constrained atoms
/// when provideds
/*!
\param ff Force field to add contributions to
\param etkdgDetails Contains information about the ETKDG force field
\param atomPairs bit set for every atom pair in the molecule where
a bit is set to one when the two atoms in the pair are distance constrained
with respect to each other
\param positions A vector of pointers to 3D Points to write out the
resulting coordinates
\param knownDistanceForceConstant force constant with which to constrain bond
distances
\param mmat Bounds matrix to use bounds from for constraints
\param numAtoms number of atoms in molecule
*/
void addLongRangeDistanceConstraints(
ForceFields::ForceField *ff,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails,
const boost::dynamic_bitset<> &atomPairs, RDGeom::Point3DPtrVect &positions,
double knownDistanceForceConstant, const BoundsMatrix &mmat,
unsigned int numAtoms) {
PRECONDITION(ff, "bad force field");
double fdist = knownDistanceForceConstant;
for (unsigned int i = 1; i < numAtoms; ++i) {
for (unsigned int j = 0; j < i; ++j) {
if (!atomPairs[j * numAtoms + i]) {
fdist = etkdgDetails.boundsMatForceScaling * 10.0;
double l = mmat.getLowerBound(i, j);
double u = mmat.getUpperBound(i, j);
if (!etkdgDetails.constrainedAtoms.empty() &&
etkdgDetails.constrainedAtoms[i] &&
etkdgDetails.constrainedAtoms[j]) {
// we're constrained, so use very tight bounds
l = u = ((*positions[i]) - (*positions[j])).length();
l -= KNOWN_DIST_TOL;
u += KNOWN_DIST_TOL;
fdist = knownDistanceForceConstant;
}
auto *contrib = new ForceFields::UFF::DistanceConstraintContrib(
ff, i, j, l, u, fdist);
ff->contribs().emplace_back(contrib);
}
}
}
}
ForceFields::ForceField *construct3DForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails) {
unsigned int N = mmat.numRows();
CHECK_INVARIANT(N == positions.size(), "");
CHECK_INVARIANT(etkdgDetails.expTorsionAtoms.size() ==
etkdgDetails.expTorsionAngles.size(),
"");
auto *field = new ForceFields::ForceField(positions[0]->dimension());
field->positions().insert(field->positions().begin(), positions.begin(),
positions.end());
// keep track which atoms are 1,2-, 1,3- or 1,4-restrained
boost::dynamic_bitset<> atomPairs(N * N);
// don't add 1-3 Distances constraints for angles where the
// central atom of the angle is the central atom of an improper torsion.
boost::dynamic_bitset<> isImproperConstrained(N);
addExperimentalTorsionTerms(field, etkdgDetails, atomPairs, N);
addImproperTorsionTerms(field, 10.0, etkdgDetails.improperAtoms,
isImproperConstrained);
add12Terms(field, etkdgDetails, atomPairs, positions,
KNOWN_DIST_FORCE_CONSTANT, N);
add13Terms(field, etkdgDetails, atomPairs, positions,
KNOWN_DIST_FORCE_CONSTANT, isImproperConstrained, true, N);
// minimum distance for all other atom pairs that aren't constrained
addLongRangeDistanceConstraints(field, etkdgDetails, atomPairs, positions,
KNOWN_DIST_FORCE_CONSTANT, mmat, N);
return field;
} // construct3DForceField
ForceFields::ForceField *construct3DForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails,
const std::map<std::pair<unsigned int, unsigned int>, double> &CPCI) {
auto *field = construct3DForceField(mmat, positions, etkdgDetails);
bool is1_4 = false;
// double dielConst = 1.0;
boost::uint8_t dielModel = 1;
for (const auto &charge : CPCI) {
auto *contrib = new ForceFields::MMFF::EleContrib(
field, charge.first.first, charge.first.second, charge.second,
dielModel, is1_4);
field->contribs().emplace_back(contrib);
}
return field;
}
ForceFields::ForceField *constructPlain3DForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails) {
unsigned int N = mmat.numRows();
CHECK_INVARIANT(N == positions.size(), "");
CHECK_INVARIANT(etkdgDetails.expTorsionAtoms.size() ==
etkdgDetails.expTorsionAngles.size(),
"");
auto *field = new ForceFields::ForceField(positions[0]->dimension());
field->positions().insert(field->positions().begin(), positions.begin(),
positions.end());
// keep track which atoms are 1,2-, 1,3- or 1,4-restrained
boost::dynamic_bitset<> atomPairs(N * N);
// don't add 1-3 Distances constraints for angles where the
// central atom of the angle is the central atom of an improper torsion.
boost::dynamic_bitset<> isImproperConstrained(N);
addExperimentalTorsionTerms(field, etkdgDetails, atomPairs, N);
add12Terms(field, etkdgDetails, atomPairs, positions,
KNOWN_DIST_FORCE_CONSTANT, N);
add13Terms(field, etkdgDetails, atomPairs, positions,
KNOWN_DIST_FORCE_CONSTANT, isImproperConstrained, false, N);
// minimum distance for all other atom pairs that aren't constrained
addLongRangeDistanceConstraints(field, etkdgDetails, atomPairs, positions,
KNOWN_DIST_FORCE_CONSTANT, mmat, N);
return field;
} // constructPlain3DForceField
ForceFields::ForceField *construct3DImproperForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const std::vector<std::vector<int>> &improperAtoms,
const std::vector<std::vector<int>> &angles,
const std::vector<int> &atomNums) {
RDUNUSED_PARAM(atomNums);
unsigned int N = mmat.numRows();
CHECK_INVARIANT(N == positions.size(), "");
auto *field = new ForceFields::ForceField(positions[0]->dimension());
field->positions().insert(field->positions().begin(), positions.begin(),
positions.end());
// improper torsions / out-of-plane bend / inversion
double oobForceScalingFactor = 10.0;
boost::dynamic_bitset<> isImproperConstrained(N);
addImproperTorsionTerms(field, oobForceScalingFactor, improperAtoms,
isImproperConstrained);
// Check that SP Centers have an angle of 180 degrees.
for (const auto &angle : angles) {
@@ -561,10 +570,9 @@ ForceFields::ForceField *construct3DImproperForceField(
auto *contrib = new ForceFields::UFF::AngleConstraintContrib(
field, angle[0], angle[1], angle[2], 179.0, 180.0,
oobForceScalingFactor);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
field->contribs().emplace_back(contrib);
}
}
return field;
} // construct3DImproperForceField
} // namespace DistGeom

29
Code/DistGeom/DistGeomUtils.h
View File

@@ -15,15 +15,13 @@
#include <Numerics/SymmMatrix.h>
#include <map>
#include <Geometry/point.h>
#include <GraphMol/ForceFieldHelpers/CrystalFF/TorsionPreferences.h>
#include "ChiralSet.h"
#include <RDGeneral/utils.h>
#include <boost/dynamic_bitset.hpp>
namespace ForceFields {
class ForceField;
namespace CrystalFF {
struct CrystalFFDetails;
}
} // namespace ForceFields
namespace DistGeom {
@@ -171,7 +169,7 @@ RDKIT_DISTGEOMETRY_EXPORT ForceFields::ForceField *constructPlain3DForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails);
//! Force field with only improper terms
//! Force field with improper terms and SP linearity contributions
/*!
\param mmat Distance bounds matrix
@@ -180,7 +178,8 @@ RDKIT_DISTGEOMETRY_EXPORT ForceFields::ForceField *constructPlain3DForceField(
\param improperAtoms A list of groups of 4 atom indices for inversion terms
\param angles List of lists with the three angle indices and whether
the center atom in the angle is SP hybridized for every angle in the molecule.
\param atomNums A list of atomic numbers for all atoms in the molecule
\param atomNums A list of atomic numbers for all atoms in the molecule,
no longer used.
\return a pointer to a ForceField with improper terms
<b>NOTE:</b> the caller is responsible for deleting this force field.
@@ -193,6 +192,26 @@ construct3DImproperForceField(
const std::vector<std::vector<int>> &angles,
const std::vector<int> &atomNums);
//! Force field with improper terms and SP linearity contributions
/*!
\param mmat Distance bounds matrix
\param positions A vector of pointers to 3D Points to write out the
resulting coordinates
\param etkdgDetails Contains information about the ETKDG force field
\return a pointer to a ForceField with improper terms
<b>NOTE:</b> the caller is responsible for deleting this force field.
*/
//! \overload
inline ForceFields::ForceField *construct3DImproperForceField(
const BoundsMatrix &mmat, RDGeom::Point3DPtrVect &positions,
const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails) {
return construct3DImproperForceField(
mmat, positions, etkdgDetails.improperAtoms, etkdgDetails.angles,
etkdgDetails.atomNums);
}
} // namespace DistGeom
#endif

5
Code/GraphMol/DistGeomHelpers/Embedder.cpp
View File

@@ -636,9 +636,8 @@ bool minimizeWithExpTorsions(RDGeom::PointPtrVect &positions,
if (embedParams.useBasicKnowledge) {
// create a force field with only the impropers
std::unique_ptr<ForceFields::ForceField> field2(
DistGeom::construct3DImproperForceField(
*eargs.mmat, positions3D, eargs.etkdgDetails->improperAtoms,
eargs.etkdgDetails->angles, eargs.etkdgDetails->atomNums));
DistGeom::construct3DImproperForceField(*eargs.mmat, positions3D,
*eargs.etkdgDetails));
if (embedParams.useRandomCoords && embedParams.coordMap != nullptr) {
for (const auto &v : *embedParams.coordMap) {
field2->fixedPoints().push_back(v.first);