[gmx-users] problem calculating PMF with distant constraints

David van der Spoel spoel at xray.bmc.uu.se
Thu Dec 15 22:04:18 CET 2005

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 
>>>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 

>>>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 
>>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 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

More information about the gromacs.org_gmx-users mailing list