From 0166c4aefc0300a40feda6669ec309fd2f610453 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Thu, 23 Apr 2026 00:21:42 -0400 Subject: [PATCH] Fix bug in inversion term for UFF, add finite difference checker. (#9228) * Fix copyright * Address review comments Removed finite diff from RDKit headers Used explicit coordinates --- Code/ForceField/CMakeLists.txt | 3 +- Code/ForceField/FiniteDifference.cpp | 53 +++++++++++++++++++++++ Code/ForceField/FiniteDifference.h | 29 +++++++++++++ Code/ForceField/UFF/CMakeLists.txt | 2 +- Code/ForceField/UFF/Inversion.cpp | 2 +- Code/ForceField/UFF/Inversions.cpp | 2 +- Code/ForceField/UFF/testUFFForceField.cpp | 40 +++++++++++++++++ 7 files changed, 127 insertions(+), 4 deletions(-) create mode 100644 Code/ForceField/FiniteDifference.cpp create mode 100644 Code/ForceField/FiniteDifference.h diff --git a/Code/ForceField/CMakeLists.txt b/Code/ForceField/CMakeLists.txt index f4f66ac81..be3e47800 100644 --- a/Code/ForceField/CMakeLists.txt +++ b/Code/ForceField/CMakeLists.txt @@ -1,6 +1,7 @@ rdkit_library(ForceField - ForceField.cpp AngleConstraint.cpp AngleConstraints.cpp + ForceField.cpp FiniteDifference.cpp + AngleConstraint.cpp AngleConstraints.cpp DistanceConstraint.cpp DistanceConstraints.cpp PositionConstraint.cpp TorsionConstraint.cpp UFF/AngleBend.cpp UFF/BondStretch.cpp UFF/Nonbonded.cpp diff --git a/Code/ForceField/FiniteDifference.cpp b/Code/ForceField/FiniteDifference.cpp new file mode 100644 index 000000000..c23143d66 --- /dev/null +++ b/Code/ForceField/FiniteDifference.cpp @@ -0,0 +1,53 @@ +// +// Copyright (C) 2026 Greg Landrum and other RDKit contributors +// +// @@ All Rights Reserved @@ +// This file is part of the RDKit. +// The contents are covered by the terms of the BSD license +// which is included in the file license.txt, found at the root +// of the RDKit source tree. +// +#include "FiniteDifference.h" +#include "ForceField.h" + +#include +#include +#include + +namespace ForceFields { + +double calcFiniteDifference(ForceField &ff, double stepSize) { + const unsigned int dim = ff.dimension(); + const unsigned int nPoints = ff.numPoints(); + const unsigned int nCoords = dim * nPoints; + + std::vector analyticGrad(nCoords, 0.0); + ff.calcGrad(analyticGrad.data()); + + std::vector pos(nCoords); + for (unsigned int i = 0; i < nPoints; ++i) { + for (unsigned int d = 0; d < dim; ++d) { + pos[i * dim + d] = (*ff.positions()[i])[d]; + } + } + + double maxDelta = 0.0; + for (unsigned int i = 0; i < nCoords; ++i) { + double orig = pos[i]; + + pos[i] = orig + stepSize; + double ePlus = ff.calcEnergy(pos.data()); + + pos[i] = orig - stepSize; + double eMinus = ff.calcEnergy(pos.data()); + + pos[i] = orig; + + double fdGrad = (ePlus - eMinus) / (2.0 * stepSize); + maxDelta = std::max(maxDelta, std::abs(fdGrad - analyticGrad[i])); + } + + return maxDelta; +} + +} // namespace ForceFields diff --git a/Code/ForceField/FiniteDifference.h b/Code/ForceField/FiniteDifference.h new file mode 100644 index 000000000..7ad04cd6a --- /dev/null +++ b/Code/ForceField/FiniteDifference.h @@ -0,0 +1,29 @@ +// +// Copyright (C) 2026 Greg Landrum and other RDKit contributors +// +// @@ All Rights Reserved @@ +// This file is part of the RDKit. +// The contents are covered by the terms of the BSD license +// which is included in the file license.txt, found at the root +// of the RDKit source tree. +// +#include +#ifndef RD_FINITEDIFFERENCE_H +#define RD_FINITEDIFFERENCE_H + +namespace ForceFields { +class ForceField; + +/// \brief Compares analytic and central-difference gradients for a force field. +/// +/// \param ff an initialized ForceField +/// \param stepSize the displacement used for the central difference +/// +/// \return the maximum absolute deviation between analytic and +/// finite-difference gradient components +RDKIT_FORCEFIELD_EXPORT double calcFiniteDifference(ForceField &ff, + double stepSize = 1e-5); + +} // namespace ForceFields + +#endif diff --git a/Code/ForceField/UFF/CMakeLists.txt b/Code/ForceField/UFF/CMakeLists.txt index f68ec4002..16befc9d5 100644 --- a/Code/ForceField/UFF/CMakeLists.txt +++ b/Code/ForceField/UFF/CMakeLists.txt @@ -1,3 +1,3 @@ rdkit_test(testUFFForceField testUFFForceField.cpp -LINK_LIBRARIES ForceFieldHelpers SmilesParse SubstructMatch +LINK_LIBRARIES ForceFieldHelpers SubstructMatch FileParsers MolTransforms ) diff --git a/Code/ForceField/UFF/Inversion.cpp b/Code/ForceField/UFF/Inversion.cpp index 729382060..0afa2ceb6 100644 --- a/Code/ForceField/UFF/Inversion.cpp +++ b/Code/ForceField/UFF/Inversion.cpp @@ -110,7 +110,7 @@ void InversionContrib::getGrad(double *pos, double *grad) const { double sinThetaSq = 1.0 - cosTheta * cosTheta; double sinTheta = std::max(sqrt(sinThetaSq), 1.0e-8); // sin(2 * W) = 2 * sin(W) * cos(W) = 2 * cos(Y) * sin(Y) - double dE_dW = -d_forceConstant * (d_C1 * cosY - 4.0 * d_C2 * cosY * sinY); + double dE_dW = -d_forceConstant * (d_C1 * cosY + 4.0 * d_C2 * cosY * sinY); RDGeom::Point3D t1 = rJL.crossProduct(rJK); RDGeom::Point3D t2 = rJI.crossProduct(rJL); RDGeom::Point3D t3 = rJK.crossProduct(rJI); diff --git a/Code/ForceField/UFF/Inversions.cpp b/Code/ForceField/UFF/Inversions.cpp index e27caec0b..d6af7d631 100644 --- a/Code/ForceField/UFF/Inversions.cpp +++ b/Code/ForceField/UFF/Inversions.cpp @@ -108,7 +108,7 @@ void InversionContribs::getGrad(double *pos, double *grad) const { 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); + (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); diff --git a/Code/ForceField/UFF/testUFFForceField.cpp b/Code/ForceField/UFF/testUFFForceField.cpp index 03424eb9a..8c2a9f415 100644 --- a/Code/ForceField/UFF/testUFFForceField.cpp +++ b/Code/ForceField/UFF/testUFFForceField.cpp @@ -14,6 +14,7 @@ #include #include +#include #include #include #include @@ -1656,6 +1657,44 @@ M END std::cerr << " done" << std::endl; } +void testFiniteDifference() { + std::cerr << "-------------------------------------" << std::endl; + std::cerr << " Test finite difference gradient check (P(F)(F)F)." + << std::endl; + + // Trigonal pyramidal PF3: P-F = 1.57 Å, F-P-F ≈ 97° + const char *molBlock = R"MOL( + RDKit 3D + + 4 3 0 0 0 0 0 0 0 0999 V2000 + 0.0000 0.0000 0.4410 P 0 0 0 0 0 0 0 0 0 0 0 0 + 1.4348 0.0000 -0.3307 F 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7174 1.2428 -0.3307 F 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7174 -1.2428 -0.3307 F 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 + 1 3 1 0 + 1 4 1 0 +M END +)MOL"; + + std::unique_ptr mol( + MolBlockToMol(molBlock, /*sanitize=*/false, /*removeHs=*/false)); + TEST_ASSERT(mol); + + std::unique_ptr ff(UFF::constructForceField(*mol)); + TEST_ASSERT(ff); + ff->initialize(); + + // Tetrahedral geometry (F-P-F ≈ 109.5°) is away from the UFF + // inversion equilibrium (w0 ≈ 84.4°), so the gradient is non-trivial. + TEST_ASSERT(ff->calcEnergy() > 1.0); + + double delta = ForceFields::calcFiniteDifference(*ff); + TEST_ASSERT(delta < 1e-6); + + std::cerr << " done" << std::endl; +} + int main() { test1(); testUFF1(); @@ -1673,4 +1712,5 @@ int main() { testUFFAllConstraints(); testUFFCopy(); testUFFButaneScan(); + testFiniteDifference(); }