cch1999 / posecheck

Pose checks for 3D Structure-based Drug Design methods
MIT License
58 stars 5 forks source link

The position constraint of 0.1 angstrom did not work as expected #8

Open AndChenCM opened 1 month ago

AndChenCM commented 1 month ago

Dear authors,

The idea to do a small local relaxation before computing strain energies to fix some inconsistency is great, but the current implementation does not seem to work as expected. You set the max displacement of atoms to be 0.1 Å using the RDKit UFFAddPositionConstraint function here, with forceConstant=1. However, if you check the displacement of atoms using the following code with a regular ligand 6ten_ligand.txt (change the suffix to sdf before testing), you will find the displacement is above the set threshold:

from rdkit import Chem
from rdkit.Chem import ChemicalForceFields

m = Chem.SDMolSupplier('6ten_ligand.sdf')[0]
m = Chem.AddHs(m, addCoords=True)
ff = ChemicalForceFields.UFFGetMoleculeForceField(m)
conf = m.GetConformer()
p = conf.GetAtomPosition(1)
ff.UFFAddPositionConstraint(1, maxDispl=0.1, forceConstant=1)
r = ff.Minimize()
q = conf.GetAtomPosition(1)
assert((p - q).Length() < 0.1)

I did some searches into this problem, and find that the force constant needs to be large enough when the displacement of the atom is larger than maxDispl, so that when the atom moves too far it gets a great energy penalty during minimization. The solution is simply changing forceConstant to 1e5 like what is done in a RDKit test script.

I have created a PR to fix this issue, along with some minor modifications. Please check if it is appropriate.