[gmx-users] Help with using g_bar

[Justin A. Lemkul]

Warren Gallin wrote:
> On 2011-04-06, at 10:58 AM, Justin A. Lemkul wrote:
> 
>>
>> Warren Gallin wrote:
>>> I am trying to use g_bar to derive a PMF curve from non-equilibrium
>>> trajectory data.  I am using v 4.5.4 on a Mac running OSX 10.6.7.
>>> I am using the end-to-end distance of a peptide as the co-ordinate of
>>> interest.  After doing a long simulation of the peptide, I selected the
>>> frames with the longest and shortest end-to-end distances, imposed a
>>> restraint on the end-to-end distance by adding a type 10 bond centered on the
>>> end-toend distance in the farme with a large force constant to the topology,
>>> and then ran a simulation to generate an ensemble of conformations with the
>>> end-to-end distances restrained at the longest and shortest distances.
>>> I then selected 5 frames from each of those two simulations and used those as
>>> the starting conformations for 20 ps simulations with free_energy set to yes
>>> with the end-to-end distance changing linearly with lambda between the two
>>> distances, with 11 equally spaced foreign_lambda values, from 0 to 1.0.  I
>>> got what I thought would be tyhe expected dhdl files (see two attached files,
>>> one for short-to-long and the other for long-to-short).
>>> However, when I invoke g_bar with these ten files as input I get no PMF, the
>>> free energy difference between the longest and shortest distances is printed
>>> out as zero, and I get a histogram file with extensive overlap.
>>> So, I must be making a mistake in either choosing the simulation parameters
>>> for the non-equilibrium trajectories or I am misunderstanding what g_bar is
>>> supposed to do, or 5 simulations in each direction is not enough to get any
>>> kind of output.  Note, I realize that I will need to run many more
>>> non-equilibrium runs to get a reliable analysis, but I am just trying to run
>>> a minimal test to be sure that I am not off on the wrong track.
>>> Could someone suggest where I am going wrong, and even better where I might
>>> find some documentation of the operational use of g_bar?
>> A complete .mdp file would help.  I'm guessing you didn't set foreign_lambda for any of these simulations?  g_bar works by evaluating energies at native and foreign lambda values.
> 
> foreign_lambda was set, hence all the ∆H curves.
> 
> I've attached the .mdp file for both directions of the simulation.
> 

OK, you've got a lot going on here:

init_lambda         =  1
delta_lambda        =  -.0001
free_energy         =  yes
foreign_lambda      =  0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

You're changing lambda as a function of time, and assigning a value of 1 to both 
foreign and native lambda.  More reliable settings have been discussed within 
the last few days, i.e.:

http://lists.gromacs.org/pipermail/gmx-users/2011-March/059694.html

So I would suggest you turn off the slow growth setting (since I am not sure how 
it will affect the BAR results, and because it's error-prone), and if that 
doesn't fix it, split the transformation into a few separate simulations at 
defined intervals of native lambda and foreign lambda (such that they don't 
overlap, either).

-Justin

> 
> 
>> A tutorial is in the works (but not for this specific application; free energy simulations in general).  No promises on when it will be done, but hopefully very soon.
>>
>> -Justin
>>
>>> Thanks,
>>> Warren Gallin
>>> ------------------------------------------------------------------------
> 
> 
> ------------------------------------------------------------------------
> 
> 
> 
> 
> ------------------------------------------------------------------------
> 
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================