Dear Rafal,
good catch: there was actually a bug in getUFFAngleBendParams(), where
the theta0 value (in radians) was multiplied by the RAD2DEG conversion
factor before being returned to the user (which is correct, since theta0
values in the UFF parameters are given in radians),but then, instead of
the original theta0 value in radians, the value in degrees was passed to
Utils::calcAngleForceConstant(), which is obviously wrong.
The bug will be fixed very soon in the GitHub repository; please find a
patch attached which will fix the problem for the time being.
Thank you very much for reporting this, kind regards
Paolo
On 04/07/2016 05:00 PM, Rafal Roszak wrote:
Dear all,
I try to reproduce UFF's Kijk (bend force constant).
In source code in the file Code/ForceField/UF/AngleBend.cpp I found
calcAngleForceConstant function. Equations in this function are the
same as in UFF Paper (JACS from 1992), however results obtained from this
function (equations) differs from this obtained using
Chem.rdForceFieldHelpers.GetUFFAngleBendParams() function in python.
Simple python script which ilustrate this is in attachment. The later
value seems to me to be more correct...
My question is: how exactly Kijk is calculated? I guess value of
calcAngleForceConstant function is modyfied somewhere (but where?)
before it is returned by Chem.rdForceFieldHelpers.GetUFFAngleBendParams()
Regards,
Rafal
------------------------------------------------------------------------------
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
diff --git a/Code/GraphMol/ForceFieldHelpers/UFF/AtomTyper.cpp
b/Code/GraphMol/ForceFieldHelpers/UFF/AtomTyper.cpp
index cf473cf..8fcf7b1 100644
--- a/Code/GraphMol/ForceFieldHelpers/UFF/AtomTyper.cpp
+++ b/Code/GraphMol/ForceFieldHelpers/UFF/AtomTyper.cpp
@@ -433,7 +433,7 @@ bool getUFFAngleBendParams(const ROMol &mol, unsigned int
idx1,
double bondOrder23 = bond[1]->getBondTypeAsDouble();
uffAngleBendParams.theta0 = RAD2DEG * paramVect[1]->theta0;
uffAngleBendParams.ka = Utils::calcAngleForceConstant(
- uffAngleBendParams.theta0, bondOrder12, bondOrder23, paramVect[0],
+ paramVect[1]->theta0, bondOrder12, bondOrder23, paramVect[0],
paramVect[1], paramVect[2]);
}
return res;
diff --git a/Code/GraphMol/ForceFieldHelpers/UFF/testUFFHelpers.cpp
b/Code/GraphMol/ForceFieldHelpers/UFF/testUFFHelpers.cpp
index f593a5f..1c62e63 100644
--- a/Code/GraphMol/ForceFieldHelpers/UFF/testUFFHelpers.cpp
+++ b/Code/GraphMol/ForceFieldHelpers/UFF/testUFFHelpers.cpp
@@ -996,7 +996,7 @@ void testUFFParamGetters() {
ForceFields::UFF::UFFAngle uffAngleBendParams;
TEST_ASSERT(UFF::getUFFAngleBendParams(*molH, 6, 7, 8,
uffAngleBendParams));
TEST_ASSERT(
- ((int)boost::math::round(uffAngleBendParams.ka * 1000) == 143325) &&
+ ((int)boost::math::round(uffAngleBendParams.ka * 1000) == 303297) &&
((int)boost::math::round(uffAngleBendParams.theta0 * 1000) == 109470));
TEST_ASSERT(
!UFF::getUFFAngleBendParams(*molH, 0, 7, 8, uffAngleBendParams));
diff --git a/Code/GraphMol/ForceFieldHelpers/Wrap/testHelpers.py
b/Code/GraphMol/ForceFieldHelpers/Wrap/testHelpers.py
index c1cea2c..7f4e0b2 100644
--- a/Code/GraphMol/ForceFieldHelpers/Wrap/testHelpers.py
+++ b/Code/GraphMol/ForceFieldHelpers/Wrap/testHelpers.py
@@ -211,7 +211,7 @@ M END"""
self.failIf(uffBondStretchParams)
uffAngleBendParams = ChemicalForceFields.GetUFFAngleBendParams(m, 6, 7, 8)
self.failUnless(uffAngleBendParams)
- self.failUnless((int(round(uffAngleBendParams[0] * 1000) == 143325))
+ self.failUnless((int(round(uffAngleBendParams[0] * 1000) == 303297))
and (int(round(uffAngleBendParams[1] * 1000) == 109470)));
uffAngleBendParams = ChemicalForceFields.GetUFFAngleBendParams(m, 0, 7, 8)
self.failIf(uffAngleBendParams)
------------------------------------------------------------------------------
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
