From a45d4d985708f6a5cb5e487581340f122773232e Mon Sep 17 00:00:00 2001 From: nmaeder <92509610+nmaeder@users.noreply.github.com> Date: Wed, 21 Aug 2024 06:06:58 +0200 Subject: [PATCH] New contribs for DG (#7711) * add angles and distances * add Inversions * add torsiona angle contribs * use new contribs in test * use new inversion and torsion contribs in dg * use new distance contribs in dg * use new angle constraints in dg * use new constraints in FF tests * update docstrings * remove unused import * include new contribs * cleanup includes * make changes requested by @greglandrum * use std::move instead of release --- Code/DistGeom/DistGeomUtils.cpp | 79 ++++++--- Code/ForceField/AngleConstraints.cpp | 143 ++++++++++++++++ Code/ForceField/AngleConstraints.h | 107 ++++++++++++ Code/ForceField/CMakeLists.txt | 13 +- Code/ForceField/DistanceConstraints.cpp | 102 ++++++++++++ Code/ForceField/DistanceConstraints.h | 103 ++++++++++++ Code/ForceField/MMFF/testMMFFForceField.cpp | 41 +++-- Code/ForceField/MMFF/testMMFFForceField.h | 1 - Code/ForceField/UFF/Contribs.h | 1 + Code/ForceField/UFF/Inversion.cpp | 66 +------- Code/ForceField/UFF/Inversion.h | 22 --- Code/ForceField/UFF/Inversions.cpp | 137 ++++++++++++++++ Code/ForceField/UFF/Inversions.h | 104 ++++++++++++ Code/ForceField/UFF/Utils.cpp | 84 ++++++++++ Code/ForceField/UFF/Utils.h | 42 +++++ Code/ForceField/UFF/testUFFForceField.cpp | 55 +++---- Code/ForceField/catch_tests.cpp | 113 +++++++++++++ .../GraphMol/ForceFieldHelpers/CMakeLists.txt | 3 +- .../CrystalFF/TorsionAngleContribs.cpp | 154 ++++++++++++++++++ .../CrystalFF/TorsionAngleContribs.h | 115 +++++++++++++ .../CrystalFF/testCrystalFF.cpp | 19 +-- .../ForceFieldHelpers/UFF/AtomTyper.cpp | 2 +- 22 files changed, 1327 insertions(+), 179 deletions(-) create mode 100644 Code/ForceField/AngleConstraints.cpp create mode 100644 Code/ForceField/AngleConstraints.h create mode 100644 Code/ForceField/DistanceConstraints.cpp create mode 100644 Code/ForceField/DistanceConstraints.h create mode 100644 Code/ForceField/UFF/Inversions.cpp create mode 100644 Code/ForceField/UFF/Inversions.h create mode 100644 Code/ForceField/UFF/Utils.cpp create mode 100644 Code/ForceField/UFF/Utils.h create mode 100644 Code/ForceField/catch_tests.cpp create mode 100644 Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionAngleContribs.cpp create mode 100644 Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionAngleContribs.h diff --git a/Code/DistGeom/DistGeomUtils.cpp b/Code/DistGeom/DistGeomUtils.cpp index 3c3a9f1a477..16579843846 100644 --- a/Code/DistGeom/DistGeomUtils.cpp +++ b/Code/DistGeom/DistGeomUtils.cpp @@ -19,11 +19,11 @@ #include #include #include -#include -#include -#include +#include +#include +#include #include -#include +#include #include #include @@ -269,6 +269,8 @@ void addImproperTorsionTerms(ForceFields::ForceField *ff, const std::vector> &improperAtoms, boost::dynamic_bitset<> &isImproperConstrained) { PRECONDITION(ff, "bad force field"); + auto inversionContribs = + std::make_unique(ff); for (const auto &improperAtom : improperAtoms) { std::vector n(4); for (unsigned int i = 0; i < 3; ++i) { @@ -293,14 +295,16 @@ void addImproperTorsionTerms(ForceFields::ForceField *ff, break; } - auto *contrib = new ForceFields::UFF::InversionContrib( - ff, improperAtom[n[0]], improperAtom[n[1]], improperAtom[n[2]], + inversionContribs->addContrib( + 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; } } + if (!inversionContribs->empty()) { + ff->contribs().push_back(std::move(inversionContribs)); + } } //! Add experimental torsion angle contributions to a force field @@ -319,6 +323,8 @@ void addExperimentalTorsionTerms( const ForceFields::CrystalFF::CrystalFFDetails &etkdgDetails, boost::dynamic_bitset<> &atomPairs, unsigned int numAtoms) { PRECONDITION(ff, "bad force field"); + auto torsionContribs = + std::make_unique(ff); for (unsigned int t = 0; t < etkdgDetails.expTorsionAtoms.size(); ++t) { int i = etkdgDetails.expTorsionAtoms[t][0]; int j = etkdgDetails.expTorsionAtoms[t][1]; @@ -329,10 +335,12 @@ void addExperimentalTorsionTerms( } 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); + torsionContribs->addContrib(i, j, k, l, + etkdgDetails.expTorsionAngles[t].second, + etkdgDetails.expTorsionAngles[t].first); + } + if (!torsionContribs->empty()) { + ff->contribs().push_back(std::move(torsionContribs)); } } @@ -355,6 +363,8 @@ void add12Terms(ForceFields::ForceField *ff, RDGeom::Point3DPtrVect &positions, double forceConstant, unsigned int numAtoms) { PRECONDITION(ff, "bad force field"); + auto distContribs = + std::make_unique(ff); for (const auto &bond : etkdgDetails.bonds) { unsigned int i = bond.first; unsigned int j = bond.second; @@ -364,9 +374,11 @@ void add12Terms(ForceFields::ForceField *ff, 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); + distContribs->addContrib(i, j, d - KNOWN_DIST_TOL, d + KNOWN_DIST_TOL, + forceConstant); + } + if (!distContribs->empty()) { + ff->contribs().push_back(std::move(distContribs)); } } //! Add 1-3 distance constraints with padding at current positions to force @@ -393,6 +405,10 @@ void add13Terms(ForceFields::ForceField *ff, const boost::dynamic_bitset<> &isImproperConstrained, bool useBasicKnowledge, unsigned int numAtoms) { PRECONDITION(ff, "bad force field"); + auto distContribs = + std::make_unique(ff); + auto angleContribs = + std::make_unique(ff); for (const auto &angle : etkdgDetails.angles) { unsigned int i = angle[0]; unsigned int j = angle[1]; @@ -405,16 +421,19 @@ void add13Terms(ForceFields::ForceField *ff, } // 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); + angleContribs->addContrib(i, j, k, 179.0, 180.0, 1); } 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); + distContribs->addContrib(i, k, d - KNOWN_DIST_TOL, d + KNOWN_DIST_TOL, + forceConstant); } } + if (!angleContribs->empty()) { + ff->contribs().push_back(std::move(angleContribs)); + } + if (!distContribs->empty()) { + ff->contribs().push_back(std::move(distContribs)); + } } //! Add long distance constraints to bounds matrix borders or constrained atoms @@ -441,6 +460,8 @@ void addLongRangeDistanceConstraints( double knownDistanceForceConstant, const BoundsMatrix &mmat, unsigned int numAtoms) { PRECONDITION(ff, "bad force field"); + auto distContribs = + std::make_unique(ff); double fdist = knownDistanceForceConstant; for (unsigned int i = 1; i < numAtoms; ++i) { for (unsigned int j = 0; j < i; ++j) { @@ -457,12 +478,13 @@ void addLongRangeDistanceConstraints( u += KNOWN_DIST_TOL; fdist = knownDistanceForceConstant; } - auto *contrib = new ForceFields::UFF::DistanceConstraintContrib( - ff, i, j, l, u, fdist); - ff->contribs().emplace_back(contrib); + distContribs->addContrib(i, j, l, u, fdist); } } } + if (!distContribs->empty()) { + ff->contribs().push_back(std::move(distContribs)); + } } ForceFields::ForceField *construct3DForceField( @@ -565,14 +587,17 @@ ForceFields::ForceField *construct3DImproperForceField( isImproperConstrained); // Check that SP Centers have an angle of 180 degrees. + auto angleContribs = + std::make_unique(field); for (const auto &angle : angles) { if (angle[3]) { - auto *contrib = new ForceFields::UFF::AngleConstraintContrib( - field, angle[0], angle[1], angle[2], 179.0, 180.0, - oobForceScalingFactor); - field->contribs().emplace_back(contrib); + angleContribs->addContrib(angle[0], angle[1], angle[2], 179.0, 180.0, + oobForceScalingFactor); } } + if (!angleContribs->empty()) { + field->contribs().push_back(std::move(angleContribs)); + } return field; } // construct3DImproperForceField } // namespace DistGeom diff --git a/Code/ForceField/AngleConstraints.cpp b/Code/ForceField/AngleConstraints.cpp new file mode 100644 index 00000000000..f519d1ce03a --- /dev/null +++ b/Code/ForceField/AngleConstraints.cpp @@ -0,0 +1,143 @@ +// +// Copyright (C) 2024 Niels Maeder 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 "AngleConstraints.h" +#include "ForceField.h" +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +namespace ForceFields { +constexpr double RAD2DEG = 180.0 / M_PI; + +AngleConstraintContribs::AngleConstraintContribs(ForceField *owner) { + PRECONDITION(owner, "bad owner"); + dp_forceField = owner; +} + +void AngleConstraintContribs::addContrib(unsigned int idx1, unsigned int idx2, + unsigned int idx3, double minAngleDeg, + double maxAngleDeg, + double forceConst) { + RANGE_CHECK(0.0, minAngleDeg, 180.0); + RANGE_CHECK(0.0, maxAngleDeg, 180.0); + URANGE_CHECK(idx1, dp_forceField->positions().size()); + URANGE_CHECK(idx2, dp_forceField->positions().size()); + URANGE_CHECK(idx3, dp_forceField->positions().size()); + PRECONDITION(maxAngleDeg >= minAngleDeg, + "minAngleDeg must be <= maxAngleDeg"); + d_contribs.emplace_back(idx1, idx2, idx3, minAngleDeg, maxAngleDeg, + forceConst); +} +void AngleConstraintContribs::addContrib(unsigned int idx1, unsigned int idx2, + unsigned int idx3, bool relative, + double minAngleDeg, double maxAngleDeg, + double forceConst) { + const RDGeom::PointPtrVect &pos = dp_forceField->positions(); + URANGE_CHECK(idx1, pos.size()); + URANGE_CHECK(idx2, pos.size()); + URANGE_CHECK(idx3, pos.size()); + PRECONDITION(maxAngleDeg >= minAngleDeg, + "minAngleDeg must be <= maxAngleDeg"); + if (relative) { + const RDGeom::Point3D &p1 = *((RDGeom::Point3D *)pos[idx1]); + const RDGeom::Point3D &p2 = *((RDGeom::Point3D *)pos[idx2]); + const RDGeom::Point3D &p3 = *((RDGeom::Point3D *)pos[idx3]); + const RDGeom::Point3D r[2] = {p1 - p2, p3 - p2}; + const double rLengthSq[2] = {std::max(1.0e-5, r[0].lengthSq()), + std::max(1.0e-5, r[1].lengthSq())}; + double cosTheta = r[0].dotProduct(r[1]) / sqrt(rLengthSq[0] * rLengthSq[1]); + cosTheta = std::clamp(cosTheta, -1.0, 1.0); + const double angle = RAD2DEG * acos(cosTheta); + minAngleDeg += angle; + maxAngleDeg += angle; + } + RANGE_CHECK(0.0, minAngleDeg, 180.0); + RANGE_CHECK(0.0, maxAngleDeg, 180.0); + d_contribs.emplace_back(idx1, idx2, idx3, minAngleDeg, maxAngleDeg, + forceConst); +} + +double AngleConstraintContribs::computeAngleTerm( + const double &angle, const AngleConstraintContribsParams &contrib) const { + double angleTerm = 0.0; + if (angle < contrib.minAngle) { + angleTerm = angle - contrib.minAngle; + } else if (angle > contrib.maxAngle) { + angleTerm = angle - contrib.maxAngle; + } + return angleTerm; +} + +double AngleConstraintContribs::getEnergy(double *pos) const { + PRECONDITION(dp_forceField, "no owner"); + PRECONDITION(pos, "bad vector"); + double accum = 0.0; + for (const auto &contrib : d_contribs) { + const RDGeom::Point3D p1(pos[3 * contrib.idx1], pos[3 * contrib.idx1 + 1], + pos[3 * contrib.idx1 + 2]); + const RDGeom::Point3D p2(pos[3 * contrib.idx2], pos[3 * contrib.idx2 + 1], + pos[3 * contrib.idx2 + 2]); + const RDGeom::Point3D p3(pos[3 * contrib.idx3], pos[3 * contrib.idx3 + 1], + pos[3 * contrib.idx3 + 2]); + const RDGeom::Point3D r[2] = {p1 - p2, p3 - p2}; + const double rLengthSq[2] = {std::max(1.0e-5, r[0].lengthSq()), + std::max(1.0e-5, r[1].lengthSq())}; + double cosTheta = r[0].dotProduct(r[1]) / sqrt(rLengthSq[0] * rLengthSq[1]); + cosTheta = std::clamp(cosTheta, -1.0, 1.0); + const double angle = RAD2DEG * acos(cosTheta); + const double angleTerm = computeAngleTerm(angle, contrib); + accum += contrib.forceConstant * angleTerm * angleTerm; + } + return accum; +} + +void AngleConstraintContribs::getGrad(double *pos, double *grad) const { + PRECONDITION(dp_forceField, "no owner"); + PRECONDITION(pos, "bad vector"); + PRECONDITION(grad, "bad vector"); + for (const auto &contrib : d_contribs) { + const RDGeom::Point3D p1(pos[3 * contrib.idx1], pos[3 * contrib.idx1 + 1], + pos[3 * contrib.idx1 + 2]); + const RDGeom::Point3D p2(pos[3 * contrib.idx2], pos[3 * contrib.idx2 + 1], + pos[3 * contrib.idx2 + 2]); + const RDGeom::Point3D p3(pos[3 * contrib.idx3], pos[3 * contrib.idx3 + 1], + pos[3 * contrib.idx3 + 2]); + const RDGeom::Point3D r[2] = {p1 - p2, p3 - p2}; + const double rLengthSq[2] = {std::max(1.0e-5, r[0].lengthSq()), + std::max(1.0e-5, r[1].lengthSq())}; + double cosTheta = r[0].dotProduct(r[1]) / sqrt(rLengthSq[0] * rLengthSq[1]); + cosTheta = std::clamp(cosTheta, -1.0, 1.0); + const double angle = RAD2DEG * acos(cosTheta); + const double angleTerm = computeAngleTerm(angle, contrib); + + const double dE_dTheta = 2.0 * RAD2DEG * contrib.forceConstant * angleTerm; + + const RDGeom::Point3D rp = r[1].crossProduct(r[0]); + const double prefactor = dE_dTheta / std::max(1.0e-5, rp.length()); + const double t[2] = {-prefactor / rLengthSq[0], prefactor / rLengthSq[1]}; + RDGeom::Point3D dedp[3]; + dedp[0] = r[0].crossProduct(rp) * t[0]; + dedp[2] = r[1].crossProduct(rp) * t[1]; + dedp[1] = -dedp[0] - dedp[2]; + double *g[3] = {&(grad[3 * contrib.idx1]), &(grad[3 * contrib.idx2]), + &(grad[3 * contrib.idx3])}; + for (unsigned int i = 0; i < 3; ++i) { + g[i][0] += dedp[i].x; + g[i][1] += dedp[i].y; + g[i][2] += dedp[i].z; + } + } +} + +} // namespace ForceFields \ No newline at end of file diff --git a/Code/ForceField/AngleConstraints.h b/Code/ForceField/AngleConstraints.h new file mode 100644 index 00000000000..0cc843e2b83 --- /dev/null +++ b/Code/ForceField/AngleConstraints.h @@ -0,0 +1,107 @@ +// +// Copyright (C) 2024 Niels Maeder 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 +#ifndef RD_ANGLECONSTRAINTS_H +#define RD_ANGLECONSTRAINTS_H +#include +#include +#include "ForceField.h" +#include "Contrib.h" + +namespace ForceFields { + +struct AngleConstraintContribsParams { + unsigned int idx1{0}; //!< index of atom1 of the angle constraint + unsigned int idx2{0}; //!< index of atom2 of the angle constraint + unsigned int idx3{0}; //!< index of atom3 of the angle constraint + double minAngle{0.0}; //!< lower bound of the flat bottom potential + double maxAngle{0.0}; //!< upper bound of the flat bottom potential + double forceConstant{1.0}; //!< force constant for angle constraint + AngleConstraintContribsParams(unsigned int idx1, unsigned int idx2, + unsigned int idx3, double minAngle, + double maxAngle, double forceConstant = 1.0) + : idx1(idx1), + idx2(idx2), + idx3(idx3), + minAngle(minAngle), + maxAngle(maxAngle), + forceConstant(forceConstant){}; +}; + +//! A term to capture all flat bottom angle constraint potentials. +class RDKIT_FORCEFIELD_EXPORT AngleConstraintContribs + : public ForceFieldContrib { + public: + AngleConstraintContribs() = default; + + //! Constructor + /*! + \param owner pointer to the owning ForceField + */ + AngleConstraintContribs(ForceField *owner); + ~AngleConstraintContribs() override = default; + //! Add a contribution to this contrib collection + /*! + \param idx1 index of atom1 in the ForceField's positions + \param idx2 index of atom2 in the ForceField's positions + \param idx3 index of atom3 in the ForceField's positions + \param minAngle minimum angle + \param maxAngle maximum angle + \param forceConst force Constant + + */ + void addContrib(unsigned int idx1, unsigned int idx2, unsigned int idx3, + double minAngleDeg, double maxAngleDeg, double forceConst); + //! Add a contribution to this contrib collection + /*! + \param idx1 index of atom1 in the ForceField's positions + \param idx2 index of atom2 in the ForceField's positions + \param idx3 index of atom3 in the ForceField's positions + \param relative whether to add the provided angle to the current angle + \param minAngle minimum angle + \param maxAngle maximum angle + \param forceConst force Constant + + */ + void addContrib(unsigned int idx1, unsigned int idx2, unsigned int idx3, + bool relative, double minAngleDeg, double maxAngleDeg, + double forceConst); + //! return the contribution of this contrib to the energy of a given state + /*! + \param pos positions of the atoms in the current state + */ + double getEnergy(double *pos) const override; + //! calculate the contribution of this contrib to the gradient at a given + /// state + /*! + \param pos positions of the atoms in the current state + \param grad gradients to be adapted + */ + void getGrad(double *pos, double *grad) const override; + //! Copy constructor + + AngleConstraintContribs *copy() const override { + return new AngleConstraintContribs(*this); + } + + //! Return true if there are no contributions in this contrib + bool empty() const { return d_contribs.empty(); } + + //! Get number of contributions in this contrib + unsigned int size() const { return d_contribs.size(); } + + private: + std::vector d_contribs; + double computeAngleTerm(const double &angle, + const AngleConstraintContribsParams &contrib) const; +}; +} // namespace ForceFields +#endif diff --git a/Code/ForceField/CMakeLists.txt b/Code/ForceField/CMakeLists.txt index a6e39114357..7c5a1f22fcc 100644 --- a/Code/ForceField/CMakeLists.txt +++ b/Code/ForceField/CMakeLists.txt @@ -1,9 +1,9 @@ -rdkit_library(ForceField +rdkit_library(ForceField AngleConstraints.cpp DistanceConstraints.cpp ForceField.cpp AngleConstraint.cpp UFF/AngleBend.cpp UFF/BondStretch.cpp UFF/Nonbonded.cpp - UFF/Inversion.cpp UFF/TorsionAngle.cpp - UFF/DistanceConstraint.cpp + UFF/Inversion.cpp UFF/Inversions.cpp UFF/TorsionAngle.cpp + UFF/DistanceConstraint.cpp UFF/Utils.cpp UFF/TorsionConstraint.cpp UFF/PositionConstraint.cpp UFF/Params.cpp MMFF/AngleBend.cpp MMFF/BondStretch.cpp MMFF/StretchBend.cpp @@ -28,7 +28,9 @@ rdkit_headers(UFF/AngleBend.h UFF/PositionConstraint.h UFF/Nonbonded.h UFF/Params.h + UFF/Utils.h UFF/Inversion.h + UFF/Inversions.h UFF/TorsionAngle.h DEST ForceField/UFF) rdkit_headers(MMFF/AngleBend.h @@ -44,6 +46,11 @@ rdkit_headers(MMFF/AngleBend.h MMFF/PositionConstraint.h MMFF/Params.h DEST ForceField/MMFF) +rdkit_catch_test(ForceFieldsCatch + catch_tests.cpp + LINK_LIBRARIES + ForceField ForceFieldHelpers) + add_subdirectory(UFF) add_subdirectory(MMFF) if(RDK_BUILD_PYTHON_WRAPPERS) diff --git a/Code/ForceField/DistanceConstraints.cpp b/Code/ForceField/DistanceConstraints.cpp new file mode 100644 index 00000000000..c5ae274c79f --- /dev/null +++ b/Code/ForceField/DistanceConstraints.cpp @@ -0,0 +1,102 @@ +// +// Copyright (C) 2024 Niels Maeder 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 "DistanceConstraints.h" +#include "ForceField.h" +#include +#include + +namespace ForceFields { +DistanceConstraintContribs::DistanceConstraintContribs(ForceField *owner) { + PRECONDITION(owner, "bad owner"); + dp_forceField = owner; +} + +void DistanceConstraintContribs::addContrib(unsigned int idx1, + unsigned int idx2, double minLen, + double maxLen, + double forceConstant) { + URANGE_CHECK(idx1, dp_forceField->positions().size()); + URANGE_CHECK(idx2, dp_forceField->positions().size()); + PRECONDITION(maxLen >= minLen, "bad bounds"); + d_contribs.emplace_back(idx1, idx2, minLen, maxLen, forceConstant); +} + +void DistanceConstraintContribs::addContrib(unsigned int idx1, + unsigned int idx2, bool relative, + double minLen, double maxLen, + double forceConstant) { + const RDGeom::PointPtrVect &pos = dp_forceField->positions(); + URANGE_CHECK(idx1, pos.size()); + URANGE_CHECK(idx2, pos.size()); + PRECONDITION(maxLen >= minLen, "bad bounds"); + if (relative) { + const RDGeom::Point3D p1 = *((RDGeom::Point3D *)pos[idx1]); + const RDGeom::Point3D p2 = *((RDGeom::Point3D *)pos[idx2]); + const auto distance = (p1 - p2).length(); + minLen = std::max(minLen + distance, 0.0); + maxLen = std::max(maxLen + distance, 0.0); + } + d_contribs.emplace_back(idx1, idx2, minLen, maxLen, forceConstant); +} + +double DistanceConstraintContribs::getEnergy(double *pos) const { + PRECONDITION(dp_forceField, "no owner"); + PRECONDITION(pos, "bad vector"); + + double accum = 0.0; + for (const auto &contrib : d_contribs) { + const auto distance2 = + dp_forceField->distance2(contrib.idx1, contrib.idx2, pos); + double difference = 0.0; + if (distance2 < contrib.minLen * contrib.minLen) { + difference = contrib.minLen - std::sqrt(distance2); + } else if (distance2 > contrib.maxLen * contrib.maxLen) { + difference = std::sqrt(distance2) - contrib.maxLen; + } else { + continue; + } + accum += 0.5 * contrib.forceConstant * difference * difference; + } + return accum; +} + +void DistanceConstraintContribs::getGrad(double *pos, double *grad) const { + PRECONDITION(dp_forceField, "no owner"); + PRECONDITION(pos, "bad vector"); + PRECONDITION(grad, "bad vector"); + + for (const auto &contrib : d_contribs) { + double preFactor = 0.0; + double distance = 0.0; + const auto distance2 = + dp_forceField->distance2(contrib.idx1, contrib.idx2, pos); + if (distance2 < contrib.minLen * contrib.minLen) { + distance = std::sqrt(distance2); + preFactor = distance - contrib.minLen; + } else if (distance2 > contrib.maxLen * contrib.maxLen) { + distance = std::sqrt(distance2); + preFactor = distance - contrib.maxLen; + } else { + continue; + } + preFactor *= contrib.forceConstant; + preFactor /= std::max(1.0e-8, distance); + const double *atom1Coords = &(pos[3 * contrib.idx1]); + const double *atom2Coords = &(pos[3 * contrib.idx2]); + for (unsigned int i = 0; i < 3; i++) { + const double dGrad = preFactor * (atom1Coords[i] - atom2Coords[i]); + grad[3 * contrib.idx1 + i] += dGrad; + grad[3 * contrib.idx2 + i] -= dGrad; + } + } +} + +} // namespace ForceFields \ No newline at end of file diff --git a/Code/ForceField/DistanceConstraints.h b/Code/ForceField/DistanceConstraints.h new file mode 100644 index 00000000000..3afa63f1d50 --- /dev/null +++ b/Code/ForceField/DistanceConstraints.h @@ -0,0 +1,103 @@ +// +// Copyright (C) 2024 Niels Maeder 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 +#ifndef RD_DISTANCECONSTRAINTS_H +#define RD_DISTANCECONSTRAINTS_H +#include +#include +#include "Contrib.h" + +namespace ForceFields { + +struct DistanceConstraintContribsParams { + unsigned int idx1{0}; //!< index of atom1 of the distance constraint + unsigned int idx2{0}; //!< index of atom2 of the distance constraint + double minLen{0.0}; //!< lower bound of the flat bottom potential + double maxLen{0.0}; //!< upper bound of the flat bottom potential + double forceConstant{1.0}; //!< force constant for distance constraint + DistanceConstraintContribsParams(unsigned int idx1, unsigned int idx2, + double minLen, double maxLen, + double forceConstant = 1.0) + : idx1(idx1), + idx2(idx2), + minLen(minLen), + maxLen(maxLen), + forceConstant(forceConstant){}; +}; + +//! A term to capture all flat bottom distance constraint potentials. +class RDKIT_FORCEFIELD_EXPORT DistanceConstraintContribs + : public ForceFieldContrib { + public: + DistanceConstraintContribs() = default; + + //! Constructor + /*! + \param owner pointer to the owning ForceField + */ + DistanceConstraintContribs(ForceField *owner); + ~DistanceConstraintContribs() override = default; + //! Add contribution to this contrib. + /*! + \param idx1 index of atom1 in the ForceField's positions + \param idx2 index of atom2 in the ForceField's positions + \param minLen minimum distance + \param maxLen maximum distance + \param forceConst force Constant + + */ + void addContrib(unsigned int idx1, unsigned int idx2, double minLen, + double maxLen, double forceConstant); + //! Add contribution to this contrib. + /*! + \param idx1 index of atom1 in the ForceField's positions + \param idx2 index of atom2 in the ForceField's positions + \param relative whether to add the provided distance to the + current distance + \param minLen minimum distance + \param maxLen maximum distance + \param forceConst force Constant + + */ + void addContrib(unsigned int idx1, unsigned int idx2, bool relative, + double minLen, double maxLen, double forceConstant); + + //! return the contribution of this contrib to the energy of a given state + /*! + \param pos positions of the atoms in the current state + */ + double getEnergy(double *pos) const override; + //! calculate the contribution of this contrib to the gradient at a given + /// state + /*! + \param pos positions of the atoms in the current state + \param grad gradients to be adapted + */ + void getGrad(double *pos, double *grad) const override; + + //! Copy constructor + DistanceConstraintContribs *copy() const override { + return new DistanceConstraintContribs(*this); + } + + //! Return true if there are no contributions in this contrib + bool empty() const { return d_contribs.empty(); } + + //! Get number of contributions in this contrib + unsigned int size() const { return d_contribs.size(); } + + private: + std::vector d_contribs; +}; + +} // namespace ForceFields + +#endif \ No newline at end of file diff --git a/Code/ForceField/MMFF/testMMFFForceField.cpp b/Code/ForceField/MMFF/testMMFFForceField.cpp index ad9289b401a..814c4a7fb1a 100644 --- a/Code/ForceField/MMFF/testMMFFForceField.cpp +++ b/Code/ForceField/MMFF/testMMFFForceField.cpp @@ -1,6 +1,5 @@ -// $Id$ // -// Copyright (C) 2013 Paolo Tosco +// Copyright (C) 2013-2024 Paolo Tosco and other RDKit contributors // // @@ All Rights Reserved @@ // This file is part of the RDKit. @@ -15,8 +14,8 @@ #include #include #include -#include -#include +#include +#include #include #include #include @@ -1449,7 +1448,7 @@ void testMMFFAllConstraints() { MMFF::MMFFMolProperties *mmffMolProperties; // distance constraints - ForceFields::MMFF::DistanceConstraintContrib *dc; + ForceFields::DistanceConstraintContribs *dc; mol = RDKit::MolBlockToMol(molBlock, true, false); TEST_ASSERT(mol); MolTransforms::setBondLength(mol->getConformer(), 1, 3, 2.0); @@ -1459,9 +1458,9 @@ void testMMFFAllConstraints() { field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties); TEST_ASSERT(field); field->initialize(); - dc = new ForceFields::MMFF::DistanceConstraintContrib(field, 1, 3, 2.0, 2.0, - 1.0e5); - field->contribs().push_back(ForceFields::ContribPtr(dc)); + dc = new ForceFields::DistanceConstraintContribs(field); + dc->addContrib(1, 3, 2.0, 2.0, 1.0e5); + field->contribs().emplace_back(dc); field->minimize(); TEST_ASSERT(RDKit::feq( MolTransforms::getBondLength(mol->getConformer(), 1, 3), 2.0, 0.1)); @@ -1472,9 +1471,9 @@ void testMMFFAllConstraints() { TEST_ASSERT(mmffMolProperties->isValid()); field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties); field->initialize(); - dc = new ForceFields::MMFF::DistanceConstraintContrib(field, 1, 3, true, -0.2, - 0.2, 1.0e5); - field->contribs().push_back(ForceFields::ContribPtr(dc)); + dc = new ForceFields::DistanceConstraintContribs(field); + dc->addContrib(1, 3, true, -0.2, 0.2, 1.0e5); + field->contribs().emplace_back(dc); field->minimize(); TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.79); delete field; @@ -1482,7 +1481,7 @@ void testMMFFAllConstraints() { delete mol; // angle constraints - ForceFields::MMFF::AngleConstraintContrib *ac; + ForceFields::AngleConstraintContribs *ac; mol = RDKit::MolBlockToMol(molBlock, true, false); TEST_ASSERT(mol); MolTransforms::setAngleDeg(mol->getConformer(), 1, 3, 6, 90.0); @@ -1492,9 +1491,9 @@ void testMMFFAllConstraints() { field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties); TEST_ASSERT(field); field->initialize(); - ac = new ForceFields::MMFF::AngleConstraintContrib(field, 1, 3, 6, 90.0, 90.0, - 100.0); - field->contribs().push_back(ForceFields::ContribPtr(ac)); + ac = new ForceFields::AngleConstraintContribs(field); + ac->addContrib(1, 3, 6, 90.0, 90.0, 100.0); + field->contribs().emplace_back(ac); field->minimize(); TEST_ASSERT(RDKit::feq( MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 90.0, 0.5)); @@ -1505,9 +1504,9 @@ void testMMFFAllConstraints() { TEST_ASSERT(mmffMolProperties->isValid()); field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties); field->initialize(); - ac = new ForceFields::MMFF::AngleConstraintContrib(field, 1, 3, 6, true, - -10.0, 10.0, 100.0); - field->contribs().push_back(ForceFields::ContribPtr(ac)); + ac = new ForceFields::AngleConstraintContribs(field); + ac->addContrib(1, 3, 6, true, -10.0, 10.0, 100.0); + field->contribs().emplace_back(ac); field->minimize(); TEST_ASSERT(RDKit::feq( MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 100.0, 0.5)); @@ -1726,9 +1725,9 @@ void testMMFFCopy() { RDKit::MMFF::constructForceField(*mol, mmffMolProperties); TEST_ASSERT(field); field->initialize(); - auto *dc = new ForceFields::MMFF::DistanceConstraintContrib( - field, 1, 3, 2.0, 2.0, 1.0e5); - field->contribs().push_back(ForceFields::ContribPtr(dc)); + auto *dc = new ForceFields::DistanceConstraintContribs(field); + dc->addContrib(1, 3, 2.0, 2.0, 1.0e5); + field->contribs().emplace_back(dc); field->minimize(); TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.99); diff --git a/Code/ForceField/MMFF/testMMFFForceField.h b/Code/ForceField/MMFF/testMMFFForceField.h index c4131c423c0..e6093d4a8ed 100644 --- a/Code/ForceField/MMFF/testMMFFForceField.h +++ b/Code/ForceField/MMFF/testMMFFForceField.h @@ -1,4 +1,3 @@ -// $Id$ // // Copyright (C) 2013 Paolo Tosco // diff --git a/Code/ForceField/UFF/Contribs.h b/Code/ForceField/UFF/Contribs.h index 1eed9681ad0..5d86693650c 100644 --- a/Code/ForceField/UFF/Contribs.h +++ b/Code/ForceField/UFF/Contribs.h @@ -17,6 +17,7 @@ #include "BondStretch.h" #include "AngleBend.h" #include "TorsionAngle.h" +#include "Inversions.h" #include "Inversion.h" #include "Nonbonded.h" diff --git a/Code/ForceField/UFF/Inversion.cpp b/Code/ForceField/UFF/Inversion.cpp index 4bb8e9f82f7..72938206033 100644 --- a/Code/ForceField/UFF/Inversion.cpp +++ b/Code/ForceField/UFF/Inversion.cpp @@ -9,6 +9,7 @@ // #include "Inversion.h" +#include "Utils.h" #include "Params.h" #include #include @@ -17,71 +18,6 @@ namespace ForceFields { namespace UFF { -namespace Utils { -double calculateCosY(const RDGeom::Point3D &iPoint, - const RDGeom::Point3D &jPoint, - const RDGeom::Point3D &kPoint, - const RDGeom::Point3D &lPoint) { - RDGeom::Point3D rJI = iPoint - jPoint; - RDGeom::Point3D rJK = kPoint - jPoint; - RDGeom::Point3D rJL = lPoint - jPoint; - rJI /= rJI.length(); - rJK /= rJK.length(); - rJL /= rJL.length(); - - RDGeom::Point3D n = rJI.crossProduct(rJK); - n /= n.length(); - - return n.dotProduct(rJL); -} - -std::tuple -calcInversionCoefficientsAndForceConstant(int at2AtomicNum, bool isCBoundToO) { - double res = 0.0; - double C0 = 0.0; - double C1 = 0.0; - double C2 = 0.0; - // if the central atom is sp2 carbon, nitrogen or oxygen - if ((at2AtomicNum == 6) || (at2AtomicNum == 7) || (at2AtomicNum == 8)) { - C0 = 1.0; - C1 = -1.0; - C2 = 0.0; - res = (isCBoundToO ? 50.0 : 6.0); - } else { - // group 5 elements are not clearly explained in the UFF paper - // the following code was inspired by MCCCS Towhee's ffuff.F - double w0 = M_PI / 180.0; - switch (at2AtomicNum) { - // if the central atom is phosphorous - case 15: - w0 *= 84.4339; - break; - - // if the central atom is arsenic - case 33: - w0 *= 86.9735; - break; - - // if the central atom is antimonium - case 51: - w0 *= 87.7047; - break; - - // if the central atom is bismuth - case 83: - w0 *= 90.0; - break; - } - C2 = 1.0; - C1 = -4.0 * cos(w0); - C0 = -(C1 * cos(w0) + C2 * cos(2.0 * w0)); - res = 22.0 / (C0 + C1 + C2); - } - res /= 3.0; - - return std::make_tuple(res, C0, C1, C2); -} -} // end of namespace Utils InversionContrib::InversionContrib(ForceField *owner, unsigned int idx1, unsigned int idx2, unsigned int idx3, diff --git a/Code/ForceField/UFF/Inversion.h b/Code/ForceField/UFF/Inversion.h index c259526d904..f895f26d4eb 100644 --- a/Code/ForceField/UFF/Inversion.h +++ b/Code/ForceField/UFF/Inversion.h @@ -52,28 +52,6 @@ class RDKIT_FORCEFIELD_EXPORT InversionContrib : public ForceFieldContrib { int d_at4Idx{-1}; double d_forceConstant, d_C0, d_C1, d_C2; }; - -namespace Utils { -//! calculates and returns the cosine of the Y angle in an improper torsion -//! (see UFF paper, equation 19) -RDKIT_FORCEFIELD_EXPORT double calculateCosY(const RDGeom::Point3D &iPoint, - const RDGeom::Point3D &jPoint, - const RDGeom::Point3D &kPoint, - const RDGeom::Point3D &lPoint); - -//! calculates and returns the UFF force constant for an improper torsion -/*! - - \param at2AtomicNum atomic number for atom 2 - \param isCBoundToO boolean flag; true if atom 2 is sp2 carbon bound to sp2 - oxygen - - \return the force constant - -*/ -RDKIT_FORCEFIELD_EXPORT std::tuple -calcInversionCoefficientsAndForceConstant(int at2AtomicNum, bool isCBoundToO); -} // namespace Utils } // namespace UFF } // namespace ForceFields #endif diff --git a/Code/ForceField/UFF/Inversions.cpp b/Code/ForceField/UFF/Inversions.cpp new file mode 100644 index 00000000000..e27caec0ba2 --- /dev/null +++ b/Code/ForceField/UFF/Inversions.cpp @@ -0,0 +1,137 @@ +// +// Copyright (C) 2024 Niels Maeder 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 "Inversions.h" +#include "Utils.h" +#include "Params.h" +#include +#include +#include +#include + +namespace ForceFields { +namespace UFF { + +InversionContribs::InversionContribs(ForceField *owner) { + PRECONDITION(owner, "bad owner"); + dp_forceField = owner; +} + +void InversionContribs::addContrib(unsigned int idx1, unsigned int idx2, + unsigned int idx3, unsigned int idx4, + int at2AtomicNum, bool isCBoundToO, + double oobForceScalingFactor) { + URANGE_CHECK(idx1, dp_forceField->positions().size()); + URANGE_CHECK(idx2, dp_forceField->positions().size()); + URANGE_CHECK(idx3, dp_forceField->positions().size()); + URANGE_CHECK(idx4, dp_forceField->positions().size()); + auto invCoeffForceCon = Utils::calcInversionCoefficientsAndForceConstant( + at2AtomicNum, isCBoundToO); + + d_contribs.emplace_back( + idx1, idx2, idx3, idx4, at2AtomicNum, isCBoundToO, + std::get<1>(invCoeffForceCon), std::get<2>(invCoeffForceCon), + std::get<3>(invCoeffForceCon), + std::get<0>(invCoeffForceCon) * oobForceScalingFactor); +} + +double InversionContribs::getEnergy(double *pos) const { + PRECONDITION(dp_forceField, "no owner"); + PRECONDITION(pos, "bad vector"); + double accum = 0; + for (const auto &contrib : d_contribs) { + const RDGeom::Point3D p1(pos[3 * contrib.idx1], pos[3 * contrib.idx1 + 1], + pos[3 * contrib.idx1 + 2]); + const RDGeom::Point3D p2(pos[3 * contrib.idx2], pos[3 * contrib.idx2 + 1], + pos[3 * contrib.idx2 + 2]); + const RDGeom::Point3D p3(pos[3 * contrib.idx3], pos[3 * contrib.idx3 + 1], + pos[3 * contrib.idx3 + 2]); + const RDGeom::Point3D p4(pos[3 * contrib.idx4], pos[3 * contrib.idx4 + 1], + pos[3 * contrib.idx4 + 2]); + const double cosY = Utils::calculateCosY(p1, p2, p3, p4); + const double sinYSq = 1.0 - cosY * cosY; + const double sinY = ((sinYSq > 0.0) ? sqrt(sinYSq) : 0.0); + // cos(2 * W) = 2 * cos(W) * cos(W) - 1 = 2 * sin(W) * sin(W) - 1 + const double cos2W = 2.0 * sinY * sinY - 1.0; + accum += contrib.forceConstant * + (contrib.C0 + contrib.C1 * sinY + contrib.C2 * cos2W); + } + return accum; +} + +void InversionContribs::getGrad(double *pos, double *grad) const { + PRECONDITION(dp_forceField, "no owner"); + PRECONDITION(pos, "bad vector"); + PRECONDITION(grad, "bad vector"); + for (const auto &contrib : d_contribs) { + const RDGeom::Point3D p1(pos[3 * contrib.idx1], pos[3 * contrib.idx1 + 1], + pos[3 * contrib.idx1 + 2]); + const RDGeom::Point3D p2(pos[3 * contrib.idx2], pos[3 * contrib.idx2 + 1], + pos[3 * contrib.idx2 + 2]); + const RDGeom::Point3D p3(pos[3 * contrib.idx3], pos[3 * contrib.idx3 + 1], + pos[3 * contrib.idx3 + 2]); + const RDGeom::Point3D p4(pos[3 * contrib.idx4], pos[3 * contrib.idx4 + 1], + pos[3 * contrib.idx4 + 2]); + double *g1 = &(grad[3 * contrib.idx1]); + double *g2 = &(grad[3 * contrib.idx2]); + double *g3 = &(grad[3 * contrib.idx3]); + double *g4 = &(grad[3 * contrib.idx4]); + RDGeom::Point3D rJI = p1 - p2; + RDGeom::Point3D rJK = p3 - p2; + RDGeom::Point3D rJL = p4 - p2; + const double dJI = rJI.length(); + const double dJK = rJK.length(); + const double dJL = rJL.length(); + if (isDoubleZero(dJI) || isDoubleZero(dJK) || isDoubleZero(dJL)) { + return; + } + rJI.normalize(); + rJK.normalize(); + rJL.normalize(); + + RDGeom::Point3D n = (-rJI).crossProduct(rJK); + n.normalize(); + double cosY = n.dotProduct(rJL); + cosY = std::clamp(cosY, -1.0, 1.0); + const double sinYSq = 1.0 - cosY * cosY; + const double sinY = std::max(sqrt(sinYSq), 1.0e-8); + double cosTheta = rJI.dotProduct(rJK); + cosTheta = std::clamp(cosTheta, -1.0, 1.0); + const double sinThetaSq = 1.0 - cosTheta * cosTheta; + const double sinTheta = std::max(sqrt(sinThetaSq), 1.0e-8); + // sin(2 * W) = 2 * sin(W) * cos(W) = 2 * cos(Y) * sin(Y) + const double dE_dW = -contrib.forceConstant * + (contrib.C1 * cosY - 4.0 * contrib.C2 * cosY * sinY); + const RDGeom::Point3D t1 = rJL.crossProduct(rJK); + const RDGeom::Point3D t2 = rJI.crossProduct(rJL); + const RDGeom::Point3D t3 = rJK.crossProduct(rJI); + const double term1 = sinY * sinTheta; + const double term2 = cosY / (sinY * sinThetaSq); + const double tg1[3] = { + (t1.x / term1 - (rJI.x - rJK.x * cosTheta) * term2) / dJI, + (t1.y / term1 - (rJI.y - rJK.y * cosTheta) * term2) / dJI, + (t1.z / term1 - (rJI.z - rJK.z * cosTheta) * term2) / dJI}; + const double tg3[3] = { + (t2.x / term1 - (rJK.x - rJI.x * cosTheta) * term2) / dJK, + (t2.y / term1 - (rJK.y - rJI.y * cosTheta) * term2) / dJK, + (t2.z / term1 - (rJK.z - rJI.z * cosTheta) * term2) / dJK}; + const double tg4[3] = {(t3.x / term1 - rJL.x * cosY / sinY) / dJL, + (t3.y / term1 - rJL.y * cosY / sinY) / dJL, + (t3.z / term1 - rJL.z * cosY / sinY) / dJL}; + for (unsigned int i = 0; i < 3; ++i) { + g1[i] += dE_dW * tg1[i]; + g2[i] += -dE_dW * (tg1[i] + tg3[i] + tg4[i]); + g3[i] += dE_dW * tg3[i]; + g4[i] += dE_dW * tg4[i]; + } + } +} +} // namespace UFF +} // namespace ForceFields \ No newline at end of file diff --git a/Code/ForceField/UFF/Inversions.h b/Code/ForceField/UFF/Inversions.h new file mode 100644 index 00000000000..db6df10115b --- /dev/null +++ b/Code/ForceField/UFF/Inversions.h @@ -0,0 +1,104 @@ +// +// Copyright (C) 2024 Niels Maeder 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 +#ifndef RD_UFFINVERSIONS_H +#define RD_UFFINVERSIONS_H +#include +#include + +namespace ForceFields { +namespace UFF { +class AtomicParams; + +struct RDKIT_FORCEFIELD_EXPORT InversionContribsParams { + unsigned int idx1{0}; //!< index of atom1 in the ForceField's positions + unsigned int idx2{0}; //!< index of atom2 in the ForceField's positions + unsigned int idx3{0}; //!< index of atom3 in the ForceField's positions + unsigned int idx4{0}; //!< index of atom4 in the ForceField's positions + int at2AtomicNum{0}; //!< atomic number for atom 2 + bool isCBoundToO{false}; //!< boolean flag; true if atom 2 is sp2 carbon + //!< bound to sp2 oxygen + double C0{0.0}; //!< inversion coefficient 0 + double C1{0.0}; //!< inversion coefficient 1 + double C2{0.0}; //!< inversion coefficient 2 + double forceConstant{1.0}; //!< force constant + InversionContribsParams(unsigned int idx1, unsigned int idx2, + unsigned int idx3, unsigned int idx4, + int at2AtomicNum, bool isCBoundToO, double C0, + double C1, double C2, double forceConstant = 1.0) + : idx1(idx1), + idx2(idx2), + idx3(idx3), + idx4(idx4), + at2AtomicNum(at2AtomicNum), + isCBoundToO(isCBoundToO), + C0(C0), + C1(C1), + C2(C2), + forceConstant(forceConstant){}; +}; +//! A term to capture all Inversion Contributionss. +class RDKIT_FORCEFIELD_EXPORT InversionContribs : public ForceFieldContrib { + public: + InversionContribs() = default; + //! Constructor + /*! + \param owner pointer to the owning ForceField + */ + InversionContribs(ForceField *owner); + + ~InversionContribs() override = default; + //! Add contribution to this contrib. + /*! + \param idx1 index of atom1 in the ForceField's positions + \param idx2 index of atom2 in the ForceField's positions + \param idx3 index of atom3 in the ForceField's positions + \param idx4 index of atom4 in the ForceField's positions + \param at2AtomicNum atomic number for atom 2 + \param isCBoundToO boolean flag; true if atom 2 is sp2 C bound to + sp2 O + \param oobForceScalingFactor scaling factor for force constant + + */ + void addContrib(unsigned int idx1, unsigned int idx2, unsigned int idx3, + unsigned int idx4, int at2AtomicNum, bool isCBoundToO, + double oobForceScalingFactor = 1.0); + //! return the contribution of this contrib to the energy of a given state + /*! + \param pos positions of the atoms in the current state + */ + double getEnergy(double *pos) const override; + //! calculate the contribution of this contrib to the gradient at a given + /// state + /*! + \param pos positions of the atoms in the current state + \param grad gradients to be adapted + */ + void getGrad(double *pos, double *grad) const override; + + //! Copy constructor + InversionContribs *copy() const override { + return new InversionContribs(*this); + } + + //! Return true if there are no contributions in this contrib + bool empty() const { return d_contribs.empty(); } + + //! Get number of contributions in this contrib + unsigned int size() const { return d_contribs.size(); } + + private: + std::vector d_contribs; +}; + +} // namespace UFF +} // namespace ForceFields + +#endif \ No newline at end of file diff --git a/Code/ForceField/UFF/Utils.cpp b/Code/ForceField/UFF/Utils.cpp new file mode 100644 index 00000000000..aff4543867e --- /dev/null +++ b/Code/ForceField/UFF/Utils.cpp @@ -0,0 +1,84 @@ +// +// Copyright (C) 2013-2024 Paolo Tosco 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 "Utils.h" +#include + +namespace ForceFields { +namespace UFF { +namespace Utils { + +double calculateCosY(const RDGeom::Point3D &iPoint, + const RDGeom::Point3D &jPoint, + const RDGeom::Point3D &kPoint, + const RDGeom::Point3D &lPoint) { + RDGeom::Point3D rJI = iPoint - jPoint; + RDGeom::Point3D rJK = kPoint - jPoint; + RDGeom::Point3D rJL = lPoint - jPoint; + rJI.normalize(); + rJK.normalize(); + rJL.normalize(); + + RDGeom::Point3D n = rJI.crossProduct(rJK); + n.normalize(); + + return n.dotProduct(rJL); +} + +std::tuple +calcInversionCoefficientsAndForceConstant(int at2AtomicNum, bool isCBoundToO) { + double res = 0.0; + double C0 = 0.0; + double C1 = 0.0; + double C2 = 0.0; + // if the central atom is sp2 carbon, nitrogen or oxygen + if ((at2AtomicNum == 6) || (at2AtomicNum == 7) || (at2AtomicNum == 8)) { + C0 = 1.0; + C1 = -1.0; + C2 = 0.0; + res = (isCBoundToO ? 50.0 : 6.0); + } else { + // group 5 elements are not clearly explained in the UFF paper + // the following code was inspired by MCCCS Towhee's ffuff.F + double w0 = M_PI / 180.0; + switch (at2AtomicNum) { + // if the central atom is phosphorous + case 15: + w0 *= 84.4339; + break; + + // if the central atom is arsenic + case 33: + w0 *= 86.9735; + break; + + // if the central atom is antimonium + case 51: + w0 *= 87.7047; + break; + + // if the central atom is bismuth + case 83: + w0 *= 90.0; + break; + } + C2 = 1.0; + C1 = -4.0 * cos(w0); + C0 = -(C1 * cos(w0) + C2 * cos(2.0 * w0)); + res = 22.0 / (C0 + C1 + C2); + } + res /= 3.0; + + return std::make_tuple(res, C0, C1, C2); +} + +} // end of namespace Utils +} // namespace UFF +} // namespace ForceFields \ No newline at end of file diff --git a/Code/ForceField/UFF/Utils.h b/Code/ForceField/UFF/Utils.h new file mode 100644 index 00000000000..bc4a1b6078d --- /dev/null +++ b/Code/ForceField/UFF/Utils.h @@ -0,0 +1,42 @@ +// +// Copyright (C) 2013-2024 Paolo Tosco 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 +#ifndef RD_UFFUTILS_H +#define RD_UFFUTILS_H +#include +#include + +namespace ForceFields { +namespace UFF { +namespace Utils { +//! calculates and returns the cosine of the Y angle in an improper torsion +//! (see UFF paper, equation 19) +RDKIT_FORCEFIELD_EXPORT double calculateCosY(const RDGeom::Point3D &iPoint, + const RDGeom::Point3D &jPoint, + const RDGeom::Point3D &kPoint, + const RDGeom::Point3D &lPoint); + +//! calculates and returns the UFF force constant for an improper torsion +/*! + + \param at2AtomicNum atomic number for atom 2 + \param isCBoundToO boolean flag; true if atom 2 is sp2 carbon bound to sp2 + oxygen + + \return the force constant + +*/ +RDKIT_FORCEFIELD_EXPORT std::tuple +calcInversionCoefficientsAndForceConstant(int at2AtomicNum, bool isCBoundToO); +} // namespace Utils +} // namespace UFF +} // namespace ForceFields + +#endif \ No newline at end of file diff --git a/Code/ForceField/UFF/testUFFForceField.cpp b/Code/ForceField/UFF/testUFFForceField.cpp index a80e0c24be2..7bc7f9a0e10 100644 --- a/Code/ForceField/UFF/testUFFForceField.cpp +++ b/Code/ForceField/UFF/testUFFForceField.cpp @@ -1,5 +1,5 @@ // -// Copyright (C) 2004-2021 Greg Landrum and other RDKit contributors +// Copyright (C) 2004-2024 Greg Landrum and other RDKit contributors // // @@ All Rights Reserved @@ // This file is part of the RDKit. @@ -21,8 +21,8 @@ #include #include #include -#include -#include +#include +#include #include #include @@ -1269,14 +1269,13 @@ void testUFFDistanceConstraints() { ff.initialize(); // C_3 - C_3, r0=1.514, k01=699.5918 - ForceFields::ForceFieldContrib *bs; - bs = new ForceFields::UFF::DistanceConstraintContrib(&ff, 0, 1, 1.35, 1.55, - 1000.0); - ff.contribs().push_back(ForceFields::ContribPtr(bs)); + auto distContribs = new ForceFields::DistanceConstraintContribs(&ff); + distContribs->addContrib(0, 1, 1.35, 1.55, 1000.0); + ff.contribs().emplace_back(distContribs); double E; - E = bs->getEnergy(p); + E = distContribs->getEnergy(p); TEST_ASSERT(RDKit::feq(E, 0.0)); - bs->getGrad(p, g); + distContribs->getGrad(p, g); for (int i = 0; i < 6; i++) { TEST_ASSERT(RDKit::feq(g[i], 0.0)); } @@ -1284,9 +1283,9 @@ void testUFFDistanceConstraints() { ff.initialize(); (*ff.positions()[1])[0] = 1.20; p[3] = 1.20; - E = bs->getEnergy(p); + E = distContribs->getEnergy(p); TEST_ASSERT(RDKit::feq(E, 11.25)); - bs->getGrad(p, g); + distContribs->getGrad(p, g); TEST_ASSERT(RDKit::feq(g[0], 150.0)); TEST_ASSERT(RDKit::feq(g[3], -150.0)); TEST_ASSERT(RDKit::feq(g[1], 0.0)); @@ -1366,50 +1365,50 @@ void testUFFAllConstraints() { ForceFields::ForceField *field; // distance constraints - ForceFields::UFF::DistanceConstraintContrib *dc; + ForceFields::DistanceConstraintContribs *dc; mol = RDKit::MolBlockToMol(molBlock, true, false); TEST_ASSERT(mol); MolTransforms::setBondLength(mol->getConformer(), 1, 3, 2.0); field = RDKit::UFF::constructForceField(*mol); TEST_ASSERT(field); field->initialize(); - dc = new ForceFields::UFF::DistanceConstraintContrib(field, 1, 3, 2.0, 2.0, - 1.0e5); - field->contribs().push_back(ForceFields::ContribPtr(dc)); + dc = new ForceFields::DistanceConstraintContribs(field); + dc->addContrib(1, 3, 2.0, 2.0, 1.0e5); + field->contribs().emplace_back(dc); field->minimize(); TEST_ASSERT(RDKit::feq( MolTransforms::getBondLength(mol->getConformer(), 1, 3), 2.0, 0.1)); delete field; field = RDKit::UFF::constructForceField(*mol); field->initialize(); - dc = new ForceFields::UFF::DistanceConstraintContrib(field, 1, 3, true, -0.2, - 0.2, 1.0e5); - field->contribs().push_back(ForceFields::ContribPtr(dc)); + dc = new ForceFields::DistanceConstraintContribs(field); + dc->addContrib(1, 3, true, -0.2, 0.2, 1.0e5); + field->contribs().emplace_back(dc); field->minimize(); TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.79); delete field; delete mol; // angle constraints - ForceFields::UFF::AngleConstraintContrib *ac; + ForceFields::AngleConstraintContribs *ac; mol = RDKit::MolBlockToMol(molBlock, true, false); TEST_ASSERT(mol); MolTransforms::setAngleDeg(mol->getConformer(), 1, 3, 6, 90.0); field = RDKit::UFF::constructForceField(*mol); TEST_ASSERT(field); field->initialize(); - ac = new ForceFields::UFF::AngleConstraintContrib(field, 1, 3, 6, 90.0, 90.0, - 100.0); - field->contribs().push_back(ForceFields::ContribPtr(ac)); + ac = new ForceFields::AngleConstraintContribs(field); + ac->addContrib(1, 3, 6, 90.0, 90.0, 100.0); + field->contribs().emplace_back(ac); field->minimize(); TEST_ASSERT(RDKit::feq( MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 90.0, 0.5)); delete field; field = RDKit::UFF::constructForceField(*mol); field->initialize(); - ac = new ForceFields::UFF::AngleConstraintContrib(field, 1, 3, 6, true, -10.0, - 10.0, 100.0); - field->contribs().push_back(ForceFields::ContribPtr(ac)); + ac = new ForceFields::AngleConstraintContribs(field); + ac->addContrib(1, 3, 6, true, -10.0, 10.0, 100.0); + field->contribs().emplace_back(ac); field->minimize(); TEST_ASSERT(RDKit::feq( MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 100.0, 0.5)); @@ -1541,9 +1540,9 @@ void testUFFCopy() { ForceFields::ForceField *field = RDKit::UFF::constructForceField(*mol); TEST_ASSERT(field); field->initialize(); - auto *dc = new ForceFields::UFF::DistanceConstraintContrib(field, 1, 3, 2.0, - 2.0, 1.0e5); - field->contribs().push_back(ForceFields::ContribPtr(dc)); + auto dc = new ForceFields::DistanceConstraintContribs(field); + dc->addContrib(1, 3, 2.0, 2.0, 1.0e5); + field->contribs().emplace_back(dc); field->minimize(); TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.99); diff --git a/Code/ForceField/catch_tests.cpp b/Code/ForceField/catch_tests.cpp new file mode 100644 index 00000000000..fd782dda749 --- /dev/null +++ b/Code/ForceField/catch_tests.cpp @@ -0,0 +1,113 @@ + +// +// Copyright (C) 2024 Niels Maeder 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 +#include +#include + +#include +#include + +#include +#include +#include + +using namespace RDKit; +namespace { + +double get_dist(const ROMol &mol, unsigned int idx1, unsigned int idx2) { + return (mol.getConformer().getAtomPos(idx1) - + mol.getConformer().getAtomPos(idx2)) + .length(); +} + +double get_angle(const ROMol &mol, unsigned int idx1, unsigned int idx2, + unsigned int idx3) { + const auto &pos1 = mol.getConformer().getAtomPos(idx1); + const auto &pos2 = mol.getConformer().getAtomPos(idx2); + const auto &pos3 = mol.getConformer().getAtomPos(idx3); + const RDGeom::Point3D r[2] = {pos1 - pos2, pos3 - pos2}; + const double rLengthSq[2] = {std::max(1.0e-5, r[0].lengthSq()), + std::max(1.0e-5, r[1].lengthSq())}; + double cosTheta = r[0].dotProduct(r[1]) / sqrt(rLengthSq[0] * rLengthSq[1]); + cosTheta = std::clamp(cosTheta, -1.0, 1.0); + return ForceFields::UFF::RAD2DEG * acos(cosTheta); +} +} // namespace + +TEST_CASE("Test DistanceConstraintContribs") { + auto mol = + "CCCO |" + "(-1.28533,-0.0567758,0.434662;-0.175447,0.695786,-0.299881;" + "0.918409,-0.342619,-0.555572;1.30936,-0.801512,0.71705)|"_smiles; + REQUIRE(mol); + SECTION("absolute distance minimization") { + auto forceField = ForceFieldsHelper::createEmptyForceFieldForMol(*mol); + REQUIRE(forceField); + forceField->initialize(); + auto contribs = std::make_unique( + forceField.get()); + contribs->addContrib(0, 1, 3.0, 3.0, 1); + contribs->addContrib(0, 2, 4.0, 4.0, 1); + forceField->contribs().push_back(std::move(contribs)); + CHECK(forceField->minimize() == 0); + CHECK(feq(get_dist(*mol, 0, 1), 3.0)); + CHECK(feq(get_dist(*mol, 0, 2), 4.0)); + } + SECTION("relative distance minimization") { + auto forceField = ForceFieldsHelper::createEmptyForceFieldForMol(*mol); + REQUIRE(forceField); + forceField->initialize(); + auto contribs = std::make_unique( + forceField.get()); + contribs->addContrib(0, 1, true, 1.0, 1.0, 1); + contribs->addContrib(0, 2, false, 4.0, 4.0, 1); + forceField->contribs().push_back(std::move(contribs)); + auto before = get_dist(*mol, 0, 1); + CHECK(forceField->minimize() == 0); + CHECK(feq(get_dist(*mol, 0, 1), 1.0 + before)); + CHECK(feq(get_dist(*mol, 0, 2), 4.0)); + } +} + +TEST_CASE("Test AngleConstraintContribs") { + auto mol = + "CCCO |" + "(-1.28533,-0.0567758,0.434662;-0.175447,0.695786,-0.299881;" + "0.918409,-0.342619,-0.555572;1.30936,-0.801512,0.71705)|"_smiles; + REQUIRE(mol); + SECTION("absolute angle minimization") { + auto forceField = ForceFieldsHelper::createEmptyForceFieldForMol(*mol); + REQUIRE(forceField); + forceField->initialize(); + auto contribs = std::make_unique( + forceField.get()); + contribs->addContrib(0, 1, 2, 120.0, 120.0, 1); + contribs->addContrib(1, 2, 3, 160.0, 160.0, 1); + forceField->contribs().push_back(std::move(contribs)); + CHECK(forceField->minimize() == 0); + CHECK(feq(get_angle(*mol, 0, 1, 2), 120.0)); + CHECK(feq(get_angle(*mol, 1, 2, 3), 160.0)); + } + SECTION("relative angle minimization") { + auto forceField = ForceFieldsHelper::createEmptyForceFieldForMol(*mol); + REQUIRE(forceField); + forceField->initialize(); + auto contribs = std::make_unique( + forceField.get()); + contribs->addContrib(0, 1, 2, true, 5.0, 5.0, 1); + contribs->addContrib(1, 2, 3, false, 160.0, 160.0, 1); + forceField->contribs().push_back(std::move(contribs)); + auto before = get_angle(*mol, 0, 1, 2); + CHECK(forceField->minimize() == 0); + CHECK(feq(get_angle(*mol, 0, 1, 2), before + 5.0)); + CHECK(feq(get_angle(*mol, 1, 2, 3), 160.0)); + } +} diff --git a/Code/GraphMol/ForceFieldHelpers/CMakeLists.txt b/Code/GraphMol/ForceFieldHelpers/CMakeLists.txt index 82b604f6c4f..1435694d40f 100644 --- a/Code/GraphMol/ForceFieldHelpers/CMakeLists.txt +++ b/Code/GraphMol/ForceFieldHelpers/CMakeLists.txt @@ -2,7 +2,7 @@ rdkit_library(ForceFieldHelpers UFF/AtomTyper.cpp UFF/Builder.cpp MMFF/AtomTyper.cpp MMFF/Builder.cpp CrystalFF/TorsionAngleM6.cpp - CrystalFF/TorsionPreferences.cpp + CrystalFF/TorsionPreferences.cpp CrystalFF/TorsionAngleContribs.cpp LINK_LIBRARIES SmilesParse SubstructMatch ForceField) target_compile_definitions(ForceFieldHelpers PRIVATE RDKIT_FORCEFIELDHELPERS_BUILD) @@ -12,6 +12,7 @@ rdkit_headers(UFF/AtomTyper.h rdkit_headers(MMFF/AtomTyper.h MMFF/Builder.h MMFF/MMFF.h DEST GraphMol/ForceFieldHelpers/MMFF) rdkit_headers(CrystalFF/TorsionAngleM6.h CrystalFF/TorsionPreferences.h + CrystalFF/TorsionAngleContribs.h DEST GraphMol/ForceFieldHelpers/CrystalFF) rdkit_catch_test(forceFieldHelpersCatch catch_tests.cpp diff --git a/Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionAngleContribs.cpp b/Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionAngleContribs.cpp new file mode 100644 index 00000000000..da94286b2f4 --- /dev/null +++ b/Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionAngleContribs.cpp @@ -0,0 +1,154 @@ +// +// Copyright (C) 2024 Niels Maeder 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 "TorsionAngleContribs.h" +#include +#include +#include +#include +#include +#include +#include + +namespace ForceFields { +namespace CrystalFF { +namespace { +double calcTorsionEnergyM6(const std::vector &forceConstants, + const std::vector &signs, const double cosPhi) { + const double cosPhi2 = cosPhi * cosPhi; + const double cosPhi3 = cosPhi * cosPhi2; + const double cosPhi4 = cosPhi * cosPhi3; + const double cosPhi5 = cosPhi * cosPhi4; + const double cosPhi6 = cosPhi * cosPhi5; + + const double cos2Phi = 2.0 * cosPhi2 - 1.0; + const double cos3Phi = 4.0 * cosPhi3 - 3.0 * cosPhi; + const double cos4Phi = 8.0 * cosPhi4 - 8.0 * cosPhi2 + 1.0; + const double cos5Phi = 16.0 * cosPhi5 - 20.0 * cosPhi3 + 5.0 * cosPhi; + const double cos6Phi = 32.0 * cosPhi6 - 48.0 * cosPhi4 + 18.0 * cosPhi2 - 1.0; + + return (forceConstants[0] * (1.0 + signs[0] * cosPhi) + + forceConstants[1] * (1.0 + signs[1] * cos2Phi) + + forceConstants[2] * (1.0 + signs[2] * cos3Phi) + + forceConstants[3] * (1.0 + signs[3] * cos4Phi) + + forceConstants[4] * (1.0 + signs[4] * cos5Phi) + + forceConstants[5] * (1.0 + signs[5] * cos6Phi)); +} +} // namespace + +TorsionAngleContribs::TorsionAngleContribs(ForceField *owner) { + PRECONDITION(owner, "bad owner"); + dp_forceField = owner; +} + +void TorsionAngleContribs::addContrib(unsigned int idx1, unsigned int idx2, + unsigned int idx3, unsigned int idx4, + std::vector forceConstants, + std::vector signs) { + PRECONDITION((idx1 != idx2) && (idx1 != idx3) && (idx1 != idx4) && + (idx2 != idx3) && (idx2 != idx4) && (idx3 != idx4), + "degenerate points"); + URANGE_CHECK(idx1, dp_forceField->positions().size()); + URANGE_CHECK(idx2, dp_forceField->positions().size()); + URANGE_CHECK(idx3, dp_forceField->positions().size()); + URANGE_CHECK(idx4, dp_forceField->positions().size()); + d_contribs.emplace_back(idx1, idx2, idx3, idx4, std::move(forceConstants), + std::move(signs)); +} + +double TorsionAngleContribs::getEnergy(double *pos) const { + PRECONDITION(dp_forceField, "no owner"); + PRECONDITION(pos, "bad vector"); + double accum = 0.0; + for (const auto &contrib : d_contribs) { + const RDGeom::Point3D iPoint(pos[3 * contrib.idx1], + pos[3 * contrib.idx1 + 1], + pos[3 * contrib.idx1 + 2]); + const RDGeom::Point3D jPoint(pos[3 * contrib.idx2], + pos[3 * contrib.idx2 + 1], + pos[3 * contrib.idx2 + 2]); + const RDGeom::Point3D kPoint(pos[3 * contrib.idx3], + pos[3 * contrib.idx3 + 1], + pos[3 * contrib.idx3 + 2]); + const RDGeom::Point3D lPoint(pos[3 * contrib.idx4], + pos[3 * contrib.idx4 + 1], + pos[3 * contrib.idx4 + 2]); + accum += calcTorsionEnergyM6( + contrib.forceConstants, contrib.signs, + MMFF::Utils::calcTorsionCosPhi(iPoint, jPoint, kPoint, lPoint)); + } + return accum; +} + +void TorsionAngleContribs::getGrad(double *pos, double *grad) const { + PRECONDITION(dp_forceField, "no owner"); + PRECONDITION(pos, "bad vector"); + PRECONDITION(grad, "bad vector"); + + for (const auto &contrib : d_contribs) { + const RDGeom::Point3D iPoint(pos[3 * contrib.idx1], + pos[3 * contrib.idx1 + 1], + pos[3 * contrib.idx1 + 2]); + const RDGeom::Point3D jPoint(pos[3 * contrib.idx2], + pos[3 * contrib.idx2 + 1], + pos[3 * contrib.idx2 + 2]); + const RDGeom::Point3D kPoint(pos[3 * contrib.idx3], + pos[3 * contrib.idx3 + 1], + pos[3 * contrib.idx3 + 2]); + const RDGeom::Point3D lPoint(pos[3 * contrib.idx4], + pos[3 * contrib.idx4 + 1], + pos[3 * contrib.idx4 + 2]); + double *g[4] = {&(grad[3 * contrib.idx1]), &(grad[3 * contrib.idx2]), + &(grad[3 * contrib.idx3]), &(grad[3 * contrib.idx4])}; + + RDGeom::Point3D r[4] = {iPoint - jPoint, kPoint - jPoint, jPoint - kPoint, + lPoint - kPoint}; + RDGeom::Point3D t[2] = {r[0].crossProduct(r[1]), r[2].crossProduct(r[3])}; + double d[2] = {t[0].length(), t[1].length()}; + if (MMFF::isDoubleZero(d[0]) || MMFF::isDoubleZero(d[1])) { + return; + } + t[0] /= d[0]; + t[1] /= d[1]; + double cosPhi = t[0].dotProduct(t[1]); + cosPhi = std::clamp(cosPhi, -1.0, 1.0); + const double sinPhiSq = 1.0 - cosPhi * cosPhi; + const double sinPhi = ((sinPhiSq > 0.0) ? sqrt(sinPhiSq) : 0.0); + const double cosPhi2 = cosPhi * cosPhi; + const double cosPhi3 = cosPhi * cosPhi2; + const double cosPhi4 = cosPhi * cosPhi3; + const double cosPhi5 = cosPhi * cosPhi4; + // dE/dPhi is independent of cartesians: + const double dE_dPhi = + (-contrib.forceConstants[0] * contrib.signs[0] * sinPhi - + 2.0 * contrib.forceConstants[1] * contrib.signs[1] * + (2.0 * cosPhi * sinPhi) - + 3.0 * contrib.forceConstants[2] * contrib.signs[2] * + (4.0 * cosPhi2 * sinPhi - sinPhi) - + 4.0 * contrib.forceConstants[3] * contrib.signs[3] * + (8.0 * cosPhi3 * sinPhi - 4.0 * cosPhi * sinPhi) - + 5.0 * contrib.forceConstants[4] * contrib.signs[4] * + (16.0 * cosPhi4 * sinPhi - 12.0 * cosPhi2 * sinPhi + sinPhi) - + 6.0 * contrib.forceConstants[4] * contrib.signs[4] * + (32.0 * cosPhi5 * sinPhi - 32.0 * cosPhi3 * sinPhi + + 6.0 * sinPhi)); + + // FIX: use a tolerance here + // this is hacky, but it's per the + // recommendation from Niketic and Rasmussen: + double sinTerm = -dE_dPhi * (MMFF::isDoubleZero(sinPhi) ? (1.0 / cosPhi) + : (1.0 / sinPhi)); + + MMFF::Utils::calcTorsionGrad(r, t, d, g, sinTerm, cosPhi); + } +} + +} // namespace CrystalFF +} // namespace ForceFields \ No newline at end of file diff --git a/Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionAngleContribs.h b/Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionAngleContribs.h new file mode 100644 index 00000000000..43673dc17c4 --- /dev/null +++ b/Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionAngleContribs.h @@ -0,0 +1,115 @@ +// +// Copyright (C) 2024 Niels Maeder 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 +#ifndef RD_TORSIONANGLECONTRIBS_H +#define RD_TORSIONANGLECONTRIBS_H +#include +#include + +namespace RDGeom { +class Point3D; +} + +namespace ForceFields { +class ForceField; +class ForceFieldContrib; +} // namespace ForceFields + +namespace ForceFields { +namespace CrystalFF { + +//! A term to capture all torsion constraint potentials. +//! +struct RDKIT_FORCEFIELDHELPERS_EXPORT TorsionAngleContribsParams { + unsigned int idx1{0}; + unsigned int idx2{0}; + unsigned int idx3{0}; + unsigned int idx4{0}; + std::vector forceConstants{6, 1.0}; + std::vector signs{6, 1}; + TorsionAngleContribsParams(unsigned int idx1, unsigned int idx2, + unsigned int idx3, unsigned int idx4, + std::vector forceConstants, + std::vector signs) + : idx1(idx1), + idx2(idx2), + idx3(idx3), + idx4(idx4), + forceConstants(forceConstants), + signs(signs) {} +}; + +class RDKIT_FORCEFIELDHELPERS_EXPORT TorsionAngleContribs + : public ForceFieldContrib { + public: + TorsionAngleContribs() = default; + + //! Constructor + /*! + \param owner pointer to the owning ForceField + */ + TorsionAngleContribs(ForceField *owner); + ~TorsionAngleContribs() = default; + //! Add contribution to this collection. + /*! + \param idx1 index of atom1 in the ForceField's positions + \param idx2 index of atom2 in the ForceField's positions + \param idx3 index of atom3 in the ForceField's positions + \param idx4 index of atom4 in the ForceField's positions + \param forceConstants force constants for the torsion potentials + \param signs phase for the torsion potentials (-1 or 1) + + */ + void addContrib(unsigned int idx1, unsigned int idx2, unsigned int idx3, + unsigned int idx4, std::vector forceConstants, + std::vector signs); + //! return the contribution of this contrib to the energy of a given state + /*! + \param pos positions of the atoms in the current state + */ + double getEnergy(double *pos) const override; + //! calculate the contribution of this contrib to the gradient at a given + /// state + /*! + \param pos positions of the atoms in the current state + \param grad gradients to be adapted + */ + void getGrad(double *pos, double *grad) const override; + //! Copy constructor + TorsionAngleContribs *copy() const override { + return new TorsionAngleContribs(*this); + } + + //! Return true if there are no contributions in this contrib + bool empty() const { return d_contribs.empty(); } + + //! Get number of contributions in this contrib + unsigned int size() const { return d_contribs.size(); } + + private: + std::vector d_contribs; +}; + +//! Calculate the torsion energy as described in 10.1021/acs.jcim.5b00654, this +//! can be used with any i > 0. +/*! + \param forceConstants Force constants for the different cosine fits + \param signs Phases of the cosine fits + \param cosPhi cosine of the torsion angle phi +*/ +RDKIT_FORCEFIELDHELPERS_EXPORT double calcTorsionEnergy( + const std::vector &forceConstants, const std::vector &signs, + const double cosPhi); + +} // namespace CrystalFF +} // namespace ForceFields + +#endif \ No newline at end of file diff --git a/Code/GraphMol/ForceFieldHelpers/CrystalFF/testCrystalFF.cpp b/Code/GraphMol/ForceFieldHelpers/CrystalFF/testCrystalFF.cpp index b1f30456bfe..2e3e2534501 100644 --- a/Code/GraphMol/ForceFieldHelpers/CrystalFF/testCrystalFF.cpp +++ b/Code/GraphMol/ForceFieldHelpers/CrystalFF/testCrystalFF.cpp @@ -15,7 +15,7 @@ #include #include #include -#include +#include #include #include #include @@ -27,7 +27,7 @@ using namespace RDGeom; using namespace RDKit; -void testTorsionAngleM6() { +void testTorsionAngleContribs() { std::cerr << "-------------------------------------" << std::endl; std::cerr << " Test CrystalFF torsional term." << std::endl; @@ -39,7 +39,6 @@ void testTorsionAngleM6() { ps.push_back(&p3); ps.push_back(&p4); - ForceFields::CrystalFF::TorsionAngleContribM6 *contrib; // ------- ------- ------- ------- ------- ------- ------- // Basic SP3 - SP3 // ------- ------- ------- ------- ------- ------- ------- @@ -49,9 +48,9 @@ void testTorsionAngleM6() { std::vector v(6, 0.0); v[2] = 4.0; - contrib = new ForceFields::CrystalFF::TorsionAngleContribM6(&ff, 0, 1, 2, 3, - v, signs); - ff.contribs().push_back(ForceFields::ContribPtr(contrib)); + auto contrib = new ForceFields::CrystalFF::TorsionAngleContribs(&ff); + contrib->addContrib(0, 1, 2, 3, v, signs); + ff.contribs().emplace_back(contrib); p1.x = 0; p1.y = 1.5; @@ -87,9 +86,9 @@ void testTorsionAngleM6() { v[1] = 7.0; ff.contribs().pop_back(); - contrib = new ForceFields::CrystalFF::TorsionAngleContribM6(&ff, 0, 1, 2, 3, - v, signs); - ff.contribs().push_back(ForceFields::ContribPtr(contrib)); + contrib = new ForceFields::CrystalFF::TorsionAngleContribs(&ff); + contrib->addContrib(0, 1, 2, 3, v, signs); + ff.contribs().emplace_back(contrib); p1.x = 0; p1.y = 1.5; @@ -221,7 +220,7 @@ int main() { BOOST_LOG(rdInfoLog) << "\t---------------------------------\n"; BOOST_LOG(rdInfoLog) << "\t SMARTS parsing\n"; - testTorsionAngleM6(); + testTorsionAngleContribs(); BOOST_LOG(rdInfoLog) << "\t---------------------------------\n"; BOOST_LOG(rdInfoLog) << "\t Seeing if non-ring torsions are applied\n"; diff --git a/Code/GraphMol/ForceFieldHelpers/UFF/AtomTyper.cpp b/Code/GraphMol/ForceFieldHelpers/UFF/AtomTyper.cpp index a0daf7c6855..97f4cece9fb 100644 --- a/Code/GraphMol/ForceFieldHelpers/UFF/AtomTyper.cpp +++ b/Code/GraphMol/ForceFieldHelpers/UFF/AtomTyper.cpp @@ -10,10 +10,10 @@ #include #include #include +#include #include #include #include -#include #include #include #include