[gmx-users] problem calculating PMF with distant constraints

[David van der Spoel]

rzangi at att.net wrote:
>  -------------- Original message ----------------------
> From: David van der Spoel <spoel at xray.bmc.uu.se>
> 
>>rzangi at att.net wrote:
>>
>>>I am trying to calculate potential of mean force between two atoms (neon, 
>>
>>methane) by running multiples simulations with distances that are constrained 
>>from lambda=0 to lambda=1 and then integrating <dHdl>. However, this doesn't 
>>seem to give the right curve as compare to calculations obtained with other 
>>methods. There appears to be a serious discrepancy.
>>
>>the dHdl curve is meaningless, only the integral is physcially 
>>meaningfull. Of course you can compare to other simulations done in an 
>>identical fashion.
> 
> 
> We compared difference free energies (i.e. the PMF curves supperimposed on each other) with identical settings and interaction parameters.
> 
There are many hidden parameters, not only in gromacs but also in other 
programs, like combination rules.
By the way you still don't say how you compute it.


> 
> 
>>>In particular when the constraint distance was bigger than the cut-off 
>>
>>interaction (so the two atoms don't `see' each other) and bigger than the range 
>>of the solvent induced interactions, the constraint force was not zero but 
>>replusive. This can be related to what was mentioned in the literature before 
>>[Chem. Phys. Lett. 196 297 (1992), JCP 100 9129(1994)] as a `dynamic' 
>>contribution to the free energy when the system allowed to rotate (the 
>>centripetal force). Thus, even when there are no forces on the atoms, but ones 
>>allow distance constraint to rotate the atoms will experience a force in the 
>>opposite direction to their center of mass, which, of course, should not be in 
>>the PMF - is there such correction to this in Gromacs?
>>
>>I don't know what to do about this, there definitely is no standard 
>>correction.
>>
>>
>>>Nevertheless, if one restrains the two atoms from moving along the x and y 
>>
>>axis and sets the constraint distance along the z-axis it still doesn't give the 
>>right answer.
>>
>>>Did someone reported these issuses before                                                                             
>>
>>?
>>I'm not completely sure that you used the right method. If you put a 
>>constraint between the atoms they won't interact at all, unless you set 
>>nrexcl to zero. An alternative method is to freeze the Neon atoms and 
>>measure the force on them.
> 
> 
> nrexcl was set to zero and fn. type 2 for the constraint `bond' was used, so they should see each other.
"should" is not good enough. have you checked the exclusions in the tpr 
file?


> 
> 
>>>This was checked with 3.2.1 and 3.1.4 versions, with a few different water
>>>models, PME and cut-offs for the electrostatic interactions as well.
>>>
>>>Is there a simple fix? The group I am working in is quite concerned about 
>>
>>this.
>>If you have a well defined system where you think you are doing exactly 
>>the same as in a publication, and you still can not reproduce it then 
>>feel free to post a bug report on bugzilla.gromacs.org. I assume you are 
>>reasonably sure that your reference values are correct.
This comment still holds. For your PMF it could e.g. be a single point.



-- 
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++