[gmx-users] problem calculating PMF with distant constraints

[rzangi at att.net]

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


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

> > 
> > 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.
> > 
> > Thanks,
> > 
> > Ronen
> > _______________________________________________
> > gmx-users mailing list
> > gmx-users at gromacs.org
> > http://www.gromacs.org/mailman/listinfo/gmx-users
> > Please don't post (un)subscribe requests to the list. Use the 
> > www interface or send it to gmx-users-request at gromacs.org.
> 
> 
> -- 
> 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
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> _______________________________________________
> gmx-users mailing list
> gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-request at gromacs.org.