From 3f3ff3858a5288fb653aa652494718648c405732 Mon Sep 17 00:00:00 2001 From: nmaeder <92509610+nmaeder@users.noreply.github.com> Date: Thu, 25 Jul 2024 06:14:12 +0200 Subject: [PATCH] Refactor distgeom minimizations (#7652) * refactor skeleton * fix contribs * adapt documentation * clarify function names * make variable names more readable and update documentation * switch boost dynamic bitset with std::vector * 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 * make changes requested by @greglandrum * make all constexprs --------- Co-authored-by: Greg Landrum --- Code/DistGeom/DistGeomUtils.cpp | 570 +++++++++++---------- Code/DistGeom/DistGeomUtils.h | 29 +- Code/GraphMol/DistGeomHelpers/Embedder.cpp | 5 +- 3 files changed, 315 insertions(+), 289 deletions(-) diff --git a/Code/DistGeom/DistGeomUtils.cpp b/Code/DistGeom/DistGeomUtils.cpp index adcb21913..3c3a9f1a4 100644 --- a/Code/DistGeom/DistGeomUtils.cpp +++ b/Code/DistGeom/DistGeomUtils.cpp @@ -28,7 +28,9 @@ #include 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 &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 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(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, 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> &improperAtoms, - const std::vector> &angles, - const std::vector &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> &improperAtoms, + boost::dynamic_bitset<> &isImproperConstrained) { + PRECONDITION(ff, "bad force field"); for (const auto &improperAtom : improperAtoms) { std::vector 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(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(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, 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> &improperAtoms, + const std::vector> &angles, + const std::vector &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 diff --git a/Code/DistGeom/DistGeomUtils.h b/Code/DistGeom/DistGeomUtils.h index ed4dfd085..65d289cf3 100644 --- a/Code/DistGeom/DistGeomUtils.h +++ b/Code/DistGeom/DistGeomUtils.h @@ -15,15 +15,13 @@ #include #include #include +#include #include "ChiralSet.h" #include #include 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 NOTE: the caller is responsible for deleting this force field. @@ -193,6 +192,26 @@ construct3DImproperForceField( const std::vector> &angles, const std::vector &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 + NOTE: 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 diff --git a/Code/GraphMol/DistGeomHelpers/Embedder.cpp b/Code/GraphMol/DistGeomHelpers/Embedder.cpp index 0bab5f8fb..eff11f7d5 100644 --- a/Code/GraphMol/DistGeomHelpers/Embedder.cpp +++ b/Code/GraphMol/DistGeomHelpers/Embedder.cpp @@ -636,9 +636,8 @@ bool minimizeWithExpTorsions(RDGeom::PointPtrVect &positions, if (embedParams.useBasicKnowledge) { // create a force field with only the impropers std::unique_ptr 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);