Add smoothing to internal-energy potentials
Explanation:
Just checked the effect of the smooth parameter on the intra energy. On a large ligand (2er7) in an fairly extended conformation (i.e. the ligand is not folding onto itself), the intra energy with smooth=0.5 is -9 kcal. However, using smooth=0.0 the intra energy is -5 kcal. This explains the different results between autodock and ocladock for ligands with many torsions.
Q:
Do you mind running some more tests on other bigger complexes too (e.g. 2er7, 3er5, 4er4) to see if that happens as well? Let me know if you can do it. Also, I assume you would need the atom-contributor pairs for each complex, wouldn't you?
By the way, I have just checked (in *.gpf and *.gpf) that I created the grids with the default smooth parameter (=0.5A) so I guess this would explain why we had no significant discrepancies on that.
A:
I don't need the pairs because I'm using the original AutoDock4.2 - I just change the smooth parameter in the .dpf file. It only affects the pairwise contributions, the grids remain unchanged.
Here are the intra energies for the complexes you asked about:
complex, smooth=0.5, smooth=0.0
2vaa, -7.64, -5.45
3er5, +12.66, +93.78
4er4, -9.84, -6.51
Q:
I see the difference. Did you implement this in autodockdevpy?. If so, I could reuse for ocladock ... thanks!
A:
Yes, it's implemented, but it's disabled in branch "ocladockenergy". Look in branch "dev".
In files pairwise_energies.py and pairwise_derivatives.py, it's implemented in function "_calc_smooth()", which modifies the distance before evaluating 'vdw' and 'hb' energy contributions.
I don't know how important it is to implement this smooth parameter. I'll try to figure it out with Stefano and gather other opinions from other people in the lab.
Q:
I have just checked the python code of "_calc_smooth()", and doesn't seem complex to implement in OpenCL. The only doubt I have is the meaning of "r" and "rij" and their relationship.
Anyway, let me know what you guys think about including this function.
A1:
"r" is a variable: it's the current distance between two atoms during the docking.
"rij" is a parameter: it's the optimum distance for the pair (e.g.: "rij" for C - C is 4 angstroms).
A2:
I just talked with Stefano and David about the smooth parameter.
The smooth parameter is important for the grids, we know that for sure. However, it is unclear how important it is for the pairwise interactions. According to the user guide (*), it was only added to pairwise interactions in version 4.2.5.
However, for the sake of publication, it would be beneficial to have a direct comparison with the current AutoDock version, so we recommend it's implementation. It would be even better if it could be a user specified argument either at run time or a compilation argument.
(*) http://autodock.scripps.edu/faqs-help/manual/autodock-4-2-user-guide/AutoDock4.2_UserGuide.pdf (see page 6)
Q:
Ok, then I will start implementing it.
A technical question: I assume "rij" (optimum distance) depends on the atoms types, doesn't it? Do you know where, either in "autodockdevpy" or AD4, I can find the "rij" values ?
A:
Yes, rij is the sum of vdW radii for the atom pair. It's calculated as 0.5 * rii + 0.5 * rjj, because rii an rjj are twice the vdW radii. The rij values must be already present in OCLADock, because they are needed to calculate C12 and C6 for vdW, and C12 and C10 for hydrogen bonds. It's probably a matter of storing them along C12, C10 and C6 for use in the energy evaluation.
Q:
Regarding the smooth parameter to be specified as a user-specified argument:
In the AD4.2 documentation, it says the force field has been optimized for a smooth value = 0.5A. So, I am setting this as the default smooth value.
But, I was wondering it such parameter has lower and upper bounds. Can you suggest these values? This would prevent any crazy smooth inputs ...
A:
in theory there wouldn't be any limits, and in the current AutoDock there isn't a check for values provided. If you want to include one, I would say that 5.0 is a pretty high upper bound, while the minimum can't go lower than 0.0 (it's a distance).