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
import math
from rdkit import Chem
from rdkit.Chem import AllChem
##all taken from UFF
rij=rjk=1.514 #natural value for C-C bond in UFF
theta=math.radians( 109.47) #natural value
zi=zk=1.912 #natural Zeff
rik = math.sqrt(rij**2 + rjk**2 - 2.0*rij*rjk*math.cos(theta) )
kijk= ((664.12*zi*zk)/(rik**5)) * (3*rij*rjk*(1- (math.cos(theta))**2) - ( rik**2 * math.cos(theta) ) )
print(" K value from equations: ")
print(kijk)
MOLEC=Chem.MolFromSmiles("CCC")
AllChem.EmbedMolecule(MOLEC)
AllChem.UFFOptimizeMolecule(MOLEC, maxIters=5540)
(k, theta)=Chem.rdForceFieldHelpers.GetUFFAngleBendParams(MOLEC, 0,1,2 )
print("\n\n K value from GetUFFAngleBendParams: ")
print(k)------------------------------------------------------------------------------
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
