[gmx-users] Re: Umbrella sampling question

[Christopher Neale]

What you reported is not what you did. It appears that grompp, gtraj, and g_dist report the same value. 
Please also report the value that you get from your pullx.xvg file that you get from mdrun, which I suspect
will also be the same.

The difference that you report is actually between the first FRAME of your trajectory from g_dist
and the first LINE of the file from the g_wham output. I see no reason to assume that the values in the
output of g_wham must be time-ordered. Also, I have never used g_wham myself (I use an external program
to do wham) and so I can not say if you are using it correctly.

My overall conclusion is that you need to investigate g_wham output not worry about a new run at this stage.

Regarding pull_pbcatom0, there is lots of information on the mailing list about this. It is a global atom number
that defines the unit cell for selection of which periodic image of each molecule will be used for the pulling.
If all of your box dimensions are >> 2*1.08 nm, then pull_pbcatom0 will not affect your results.

Chris.

-- original message --

Gmx QA gmxquestions at gmail.com 
Wed Nov 14 15:06:46 CET 2012
Previous message: [gmx-users] Weird result of WHAM
Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
Hi Chris,

and thank you for your reply. I should have included my g_dist command in
my first mail. Here goes:

I first run trjconv to extract the individual frames from my pulling
trajectory
$ trjconv -f pull.xtc -s pull.tpr -o conf.gro -sep -pbc mol -ur compact

Then, g_dist like so:
$ g_dist -s pull.tpr -f conf0.gro -n index.ndx -o dist0.xvg

where index.ndx is an index-file that contains both pull-groups (ie the
Protein in one group, and another group with the single atom from the
molecule I'm pulling). I've checked that is is the same as in the mdp-file
for the pulling simulation.

>From this dist0.xvg-file I extract the z-component of the distance:
$ tail -n 1 dist0.xvg | awk '{print $5}'
-1.8083301

And this is also the same value I get using g_traj to manually calculate it.

Now, I run g_wham with this command line:
$ g_wham -it tpr-files.dat -if pullf-files.dat -o -hist -unit kCal

which gives me a profile.xvg file (per the default -o flag) containing the
PMF.
The top-part of the profile.xvg file looks like this:

# This file was created Tue Nov  6 09:56:23 2012
# by the following command:
# g_wham -it tpr-files.dat -if pullf-files.dat -o -hist -unit kCal
#
# g_wham is part of G R O M A C S:
#
# Gromacs Runs On Most of All Computer Systems
#
@    title "Umbrella potential"
@    xaxis  label "z"
@    yaxis  label "E (kcal mol\S-1\N)"
@TYPE xy
1.761971e+00    0.000000e+00
1.782558e+00    -1.436057e-01
1.803145e+00    -3.857955e-01

and I assumed that the distances reported here should be the same as the
distances calculated by g_dist for each frame, but as you can see they are
not (1.761971 vs 1.8083301)

I have been trying to understand how the 1.761971 distance is calculated,
but I do not get it since I use the pullf.xvg-files as input to g_wham. And
even using the z-component of the distance from the pullx.xvg-file and
calculating ie the average does not give me a value close to 1.761971. For
the next couple of frames, there are also differences, but they are not
constant.


I should add that after having sent my first mail, I found this thread:
http://gromacs.5086.n6.nabble.com/umbrella-sampling-PMF-position-discrepancy-td5000886.html
which fairly well describes my problem. Following the suggestion to change
(or in my case add, it was not there originally) the pull_pbcatom0 entry to
the mdp-file as an atom that is close to the COM of the protein, gives me a
distance of 1.809 reported from running grompp (when making the tpr for the
pulling simulation). I've started another pulling run with this new tpr,
but I don't like making changes I don't understand. I guess I'm not clear
about what pull_pbcatom0 does, but I though the distance in this case was
rather unambiguous.

Thanks
/PK