Torsion constraints not working properly for negative torsions

[greglandrum]

Short summary: torsion constraints do not seem to work if the dihedral either starts negative or becomes negative.

In [22]:

m = Chem.MolFromSmiles('CCCC')
AllChem.EmbedMolecule(m,randomSeed=0xf00d)
confid=0
conf = cm.GetConformer(confid)

# we can't just create an empty force field from python, so build a version of the
# molecule that has no bonds:
AllChem.SetDihedralDeg(m.GetConformer(confid),0,1,2,3,50)
tm = Chem.EditableMol(m)
for i in range(m.GetNumBonds()):
    bnd = m.GetBondWithIdx(i)
    tm.RemoveBond(bnd.GetBeginAtomIdx(),bnd.GetEndAtomIdx())
nm = tm.GetMol()

# build a force field for that no-bond molecule by constraining all distances
# that aren't the torsion we care about to their current value += .01
# in a real molecule, we'd probably want to only do 1-2 and 1-3 distances 
# (i.e. bonds and angles) here in order to avoid over-constraining the system
dm = AllChem.Get3DDistanceMatrix(nm,confId=confid)


ff = rdff.UFFGetMoleculeForceField(nm,confId=confid)
ff.AddDistanceConstraint(0,1,dm[0,1]-.01,dm[0,1]+.01,1000)
ff.AddDistanceConstraint(0,2,dm[0,2]-.01,dm[0,2]+.01,1000)
ff.AddDistanceConstraint(1,2,dm[1,2]-.01,dm[1,2]+.01,1000)
ff.AddDistanceConstraint(1,3,dm[1,3]-.01,dm[1,3]+.01,1000)
ff.AddDistanceConstraint(2,3,dm[2,3]-.01,dm[2,3]+.01,1000)

# add the torsion contraint:
ff.UFFAddTorsionConstraint(0,1,2,3,False,10,20,1000)
print('starting positive')
print("init: ",ff.CalcEnergy())
print("   dihed1: ",AllChem.GetDihedralDeg(nm.GetConformer(confid),0,1,2,3))

ff.Minimize()

print("done: ",ff.CalcEnergy())
print("   dihed1: ",AllChem.GetDihedralDeg(nm.GetConformer(confid),0,1,2,3))


# now repeat that, but set the initial dihedral negative:

AllChem.SetDihedralDeg(m.GetConformer(confid),0,1,2,3,-30)
tm = Chem.EditableMol(m)
for i in range(m.GetNumBonds()):
    bnd = m.GetBondWithIdx(i)
    tm.RemoveBond(bnd.GetBeginAtomIdx(),bnd.GetEndAtomIdx())
nm = tm.GetMol()


ff = rdff.UFFGetMoleculeForceField(nm,confId=confid)
ff.AddDistanceConstraint(0,1,dm[0,1]-.01,dm[0,1]+.01,1000)
ff.AddDistanceConstraint(0,2,dm[0,2]-.01,dm[0,2]+.01,1000)
ff.AddDistanceConstraint(1,2,dm[1,2]-.01,dm[1,2]+.01,1000)
ff.AddDistanceConstraint(1,3,dm[1,3]-.01,dm[1,3]+.01,1000)
ff.AddDistanceConstraint(2,3,dm[2,3]-.01,dm[2,3]+.01,1000)
# add the torsion contraint:
ff.UFFAddTorsionConstraint(0,1,2,3,False,-10,-8,1000)
print('starting negative')
print("init: ",ff.CalcEnergy())
print("   dihed1: ",AllChem.GetDihedralDeg(nm.GetConformer(confid),0,1,2,3))

ff.Minimize()

print("done: ",ff.CalcEnergy())
print("   dihed1: ",AllChem.GetDihedralDeg(nm.GetConformer(confid),0,1,2,3))

# initial dihedral positive, constraint negative:

AllChem.SetDihedralDeg(m.GetConformer(confid),0,1,2,3,30)
tm = Chem.EditableMol(m)
for i in range(m.GetNumBonds()):
    bnd = m.GetBondWithIdx(i)
    tm.RemoveBond(bnd.GetBeginAtomIdx(),bnd.GetEndAtomIdx())
nm = tm.GetMol()


ff = rdff.UFFGetMoleculeForceField(nm,confId=confid)
ff.AddDistanceConstraint(0,1,dm[0,1]-.01,dm[0,1]+.01,1000)
ff.AddDistanceConstraint(0,2,dm[0,2]-.01,dm[0,2]+.01,1000)
ff.AddDistanceConstraint(1,2,dm[1,2]-.01,dm[1,2]+.01,1000)
ff.AddDistanceConstraint(1,3,dm[1,3]-.01,dm[1,3]+.01,1000)
ff.AddDistanceConstraint(2,3,dm[2,3]-.01,dm[2,3]+.01,1000)
# add the torsion contraint:
ff.UFFAddTorsionConstraint(0,1,2,3,False,-10,-8,1000)
print('pos -> neg')
print("init: ",ff.CalcEnergy())
print("   dihed1: ",AllChem.GetDihedralDeg(nm.GetConformer(confid),0,1,2,3))

ff.Minimize()

print("done: ",ff.CalcEnergy())
print("   dihed1: ",AllChem.GetDihedralDeg(nm.GetConformer(confid),0,1,2,3))

produces:

starting positive
init:  137.07783890401893
   dihed1:  50.00000000000001
done:  1.1286835858560129e-17
   dihed1:  19.597109864499913
starting negative
init:  60.923483957341716
   dihed1:  -29.999999999999996
done:  60.923483957341716
   dihed1:  -29.999999999999996
pos -> neg
init:  219.9337770860036
   dihed1:  30.00000000000001
done:  49.00453730195457
   dihed1:  -23.019269622710254

[greglandrum]

Fixed as part of PR #395