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

@@ -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

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 <algorithm>
+#include <cmath>
+#include <vector>
+
+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<double> analyticGrad(nCoords, 0.0);
+  ff.calcGrad(analyticGrad.data());
+
+  std::vector<double> 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

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 <RDGeneral/export.h>
+#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

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 )

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);

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);

Code/ForceField/UFF/testUFFForceField.cpp

@@ -14,6 +14,7 @@
 #include <RDGeneral/utils.h>
 #include <Geometry/point.h>
 
+#include <ForceField/FiniteDifference.h>
 #include <ForceField/ForceField.h>
 #include <ForceField/UFF/Params.h>
 #include <ForceField/UFF/BondStretch.h>
@@ -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<ROMol> mol(
+      MolBlockToMol(molBlock, /*sanitize=*/false, /*removeHs=*/false));
+  TEST_ASSERT(mol);
+
+  std::unique_ptr<ForceFields::ForceField> 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();
 }