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

Fixes Github #7901 (#8264)

Browse Source 
* does not work

* basically works
(at least it seems to)
needs more tesing and verification

* optimization and improved tests

passes local tests

* response to review
This commit is contained in:
Greg Landrum
2025-02-13 12:46:56 +01:00
committed by greg landrum
parent c4ca101094
commit 497c32a147
4 changed files with 111 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

51
Code/ForceField/UFF/AngleBend.cpp
View File

@@ -19,6 +19,10 @@
namespace ForceFields {
namespace UFF {
namespace {
constexpr double ANGLE_CORRECTION_THRESHOLD = 0.8660;
}
namespace Utils {
double calcAngleForceConstant(double theta0, double bondOrder12,
double bondOrder23, const AtomicParams *at1Params,
@@ -81,20 +85,20 @@ AngleBendContrib::AngleBendContrib(ForceField *owner, unsigned int idx1,
URANGE_CHECK(idx3, owner->positions().size());
// the following is a hack to get decent geometries
// with 3- and 4-membered rings incorporating sp2 atoms
double theta0 = at2Params->theta0;
d_theta0 = at2Params->theta0;
if (order >= 30) {
switch (order) {
case 30:
theta0 = 150.0 / 180.0 * M_PI;
d_theta0 = 150.0 / 180.0 * M_PI;
break;
case 35:
theta0 = 60.0 / 180.0 * M_PI;
d_theta0 = 60.0 / 180.0 * M_PI;
break;
case 40:
theta0 = 135.0 / 180.0 * M_PI;
d_theta0 = 135.0 / 180.0 * M_PI;
break;
case 45:
theta0 = 90.0 / 180.0 * M_PI;
d_theta0 = 90.0 / 180.0 * M_PI;
break;
}
order = 0;
@@ -106,10 +110,10 @@ AngleBendContrib::AngleBendContrib(ForceField *owner, unsigned int idx1,
d_at3Idx = idx3;
d_order = order;
d_forceConstant = Utils::calcAngleForceConstant(
theta0, bondOrder12, bondOrder23, at1Params, at2Params, at3Params);
d_theta0, bondOrder12, bondOrder23, at1Params, at2Params, at3Params);
if (order == 0) {
double sinTheta0 = sin(theta0);
double cosTheta0 = cos(theta0);
double sinTheta0 = sin(d_theta0);
double cosTheta0 = cos(d_theta0);
d_C2 = 1. / (4. * std::max(sinTheta0 * sinTheta0, 1e-8));
d_C1 = -4. * d_C2 * cosTheta0;
d_C0 = d_C2 * (2. * cosTheta0 * cosTheta0 + 1.);
@@ -139,6 +143,19 @@ double AngleBendContrib::getEnergy(double *pos) const {
double angleTerm = getEnergyTerm(cosTheta, sinThetaSq);
double res = d_forceConstant * angleTerm;
// The original UFF does not include any penalty for angles that are zero
// degrees.
// This can lead to overlapping 1-3 atoms (e.g.e Github #7901), which is
// obviously bad. We add an empiricial penalty for angles close to zero
// borrowed from OpenBabel such that the energy goes up exponentially if the
// angle is less than approx theta0,
// For the sake of efficiency, we only add the penalty if the angle is less
// than 30 degrees
if (d_order && d_order < 5 && cosTheta > ANGLE_CORRECTION_THRESHOLD) {
auto theta = acos(cosTheta);
res += exp(-20.0 * (theta - d_theta0 + 0.25));
}
return res;
}
@@ -164,15 +181,27 @@ void AngleBendContrib::getGrad(double *pos, double *grad) const {
double sinThetaSq = 1.0 - cosTheta * cosTheta;
double sinTheta = std::max(sqrt(sinThetaSq), 1.0e-8);
// std::cerr << "GRAD: " << cosTheta << " (" << acos(cosTheta)<< "), ";
// std::cerr << sinTheta << " (" << asin(sinTheta)<< ")" << std::endl;
// use the chain rule:
// dE/dx = dE/dTheta * dTheta/dx
// dE/dTheta is independent of cartesians:
double dE_dTheta = getThetaDeriv(cosTheta, sinTheta);
// The original UFF does not include any penalty for angles that are zero
// degrees.
// This can lead to overlapping 1-3 atoms (e.g.e Github #7901), which is
// obviously bad. We add an empiricial penalty for angles close to zero
// borrowed from OpenBabel such that the energy goes up exponentially if
// the angle is less than approx theta0
// For the sake of efficiency, we only add the penalty if the angle is less
// than 30 degrees
if (d_order && d_order < 5 && cosTheta > ANGLE_CORRECTION_THRESHOLD) {
auto theta = acos(cosTheta);
auto corr = -20.0 * exp(-20.0 * (theta - d_theta0 + 0.25));
dE_dTheta += corr;
}
Utils::calcAngleBendGrad(r, dist, g, dE_dTheta, cosTheta, sinTheta);
}

4
Code/ForceField/UFF/AngleBend.h
View File

@@ -59,8 +59,8 @@ class RDKIT_FORCEFIELD_EXPORT AngleBendContrib : public ForceFieldContrib {
int d_at2Idx{-1};
int d_at3Idx{-1};
unsigned int d_order{0};
double d_forceConstant, d_C0, d_C1, d_C2;
double d_forceConstant, d_C0, d_C1, d_C2, d_theta0;
double getEnergyTerm(double cosTheta, double sinThetaSq) const;
double getThetaDeriv(double cosTheta, double sinTheta) const;
};

2
Code/GraphMol/ForceFieldHelpers/CMakeLists.txt
View File

@@ -16,7 +16,7 @@ rdkit_headers(CrystalFF/TorsionAngleM6.h CrystalFF/TorsionPreferences.h
DEST GraphMol/ForceFieldHelpers/CrystalFF)
rdkit_catch_test(forceFieldHelpersCatch catch_tests.cpp
LINK_LIBRARIES ForceFieldHelpers SmilesParse ForceField )
LINK_LIBRARIES MolTransforms ForceFieldHelpers SmilesParse FileParsers ForceField )
add_subdirectory(MMFF)
add_subdirectory(UFF)

68
Code/GraphMol/ForceFieldHelpers/catch_tests.cpp
View File

@@ -13,8 +13,11 @@
#include <catch2/catch_all.hpp>
#include <GraphMol/SmilesParse/SmilesParse.h>
#include <GraphMol/FileParsers/FileParsers.h>
#include <GraphMol/ForceFieldHelpers/UFF/UFF.h>
#include <ForceField/MMFF/Params.h>
#include <ForceField/MMFF/BondStretch.h>
#include <GraphMol/MolTransforms/MolTransforms.h>
#include "FFConvenience.h"
@@ -52,3 +55,68 @@ TEST_CASE("Test empty force field") {
CHECK(std::round(dist) == 100);
}
}
TEST_CASE("github #7901") {
#if 1
auto mb = R"CTAB(Acetone, enolate form
RDKit 3D
9 8 0 0 0 0 0 0 0 0999 V2000
-1.2143 -0.3494 0.0962 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1072 0.3398 0.0041 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1682 1.7287 0.0554 O 0 0 0 0 0 1 0 0 0 0 0 0
1.2210 -0.3737 -0.1263 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.2311 -1.0203 0.9811 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.3951 -0.9468 -0.8225 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.0281 0.3987 0.2006 H 0 0 0 0 0 0 0 0 0 0 0 0
2.1862 0.1115 -0.1943 H 0 0 0 0 0 0 0 0 0 0 0 0
2.8861 0.1115 -0.1943 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
2 4 2 0
1 5 1 0
1 6 1 0
1 7 1 0
4 8 1 0
4 9 1 0
M CHG 1 3 -1
M END)CTAB";
#else
auto mb = R"CTAB(Acetone, enolate form
RDKit 3D
9 8 0 0 0 0 0 0 0 0999 V2000
-1.2143 -0.3494 0.0962 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1072 0.3398 0.0041 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1682 1.7287 0.0554 O 0 0 0 0 0 1 0 0 0 0 0 0
1.2210 -0.3737 -0.1263 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.2311 -1.0203 0.9811 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.3951 -0.9468 -0.8225 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.0281 0.3987 0.2006 H 0 0 0 0 0 0 0 0 0 0 0 0
2.1862 0.1115 -0.1943 H 0 0 0 0 0 0 0 0 0 0 0 0
2.8861 -1.3115 -0.1943 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
2 4 2 0
1 5 1 0
1 6 1 0
1 7 1 0
4 8 1 0
4 9 1 0
M CHG 1 3 -1
M END)CTAB";
#endif
v2::FileParsers::MolFileParserParams params;
params.removeHs = false;
auto mol = v2::FileParsers::MolFromMolBlock(mb, params);
REQUIRE(mol);
auto &conf = mol->getConformer();
std::unique_ptr<ForceFields::ForceField> ff{UFF::constructForceField(*mol)};
ff->initialize();
auto needsMore = ff->minimize(200);
CHECK(!needsMore);
CHECK((conf.getAtomPos(7) - conf.getAtomPos(8)).length() > 1.80);
CHECK(MolTransforms::getAngleDeg(conf, 1, 3, 8) > 110);
CHECK(MolTransforms::getAngleDeg(conf, 1, 3, 7) > 110);
CHECK(MolTransforms::getAngleDeg(conf, 7, 3, 8) > 110);
}