[gmx-users] Re: Umbrella sampling question

Wed Nov 14 15:06:46 CET 2012

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

Can you please let us know exactly how you got the two values that you find
> to be different (but expected
> to be the same)? i.e. post your full g_dist command and explain how you
> observed the value in the output from
> mdrun. One frame should be enough for now (as long as you are sure -- and
> can prove to us -- that it is
> indeed the same frame).
>
> You don't define profile.xvg and I wonder if the discrepancy might be
> related to your use of pull_start=yes
> Do you see a consistent bias between mdrun output and g_dist output if you
> subtract the values for
> different time frames? If so, is this difference the same as the distance
> in the initial frame (pull_start=yes)?
>
> Chris.
>
> -- original message --
>
> I'm a new poster on the maillist, and new to umbrella sampling but not to
> MD in general.
>
> I have recently done some pulling simulation with Gromacs, but have a few
> questions about the outcome, and in particular the distances that are
> calculated for the different windows. I realize this is a question that has
> come up before, and there are some useful posts in the archives, but
> nothing that exactly answers my questions.
>
> To begin, I ran a pulling-simulation pulling a molecule starting in the
> middle of a membrane (and inside a protein transport channel) straight up
> in z. To that end, I used the following mdp-settings:
>
> pull = umbrella
> pull_geometry = direction
> pull_vec1 = 0 0 1
> pull_ngroups = 1
> pull_start = yes
> pull_group0 = Protein
> pull_group1 = pg1
> pull_rate1 = 0.001
> pull_k1 = 500
>
> So my pull groups are the COM of the protein, and one atom named pg1 on the
> pull molecule.
> To my understanding, this should pull the pg1-molecule straight up in z,
> and that is indeed also what is happening in the simulation.
>
> I ran this simulation for 6 ns, and it resulted in about 40 separate
> conformations to use for umbrella sampling. All of those simulations also
> seem to work, I can use the resulting pullf,xvg-files as input to g_wham,
> and get a histogram-plot with reasonable overlaps.
>
> However, I'm trying to understand how the various distances relate to each
> other. For example, in the profile.xvg-file I get a z-distance for the
> first frame as 1.761971 nm, while checking with g_dist gives me 1.80833
> (for z). Continuing a few frames further, there are still differences, and
> they appear to be random.
>
> How are the distances in the profile.xvg-file computed? The average of the
> dZ column of the pullx-file for the first frame is 1.75787, which is sort
> of close to 1.761971, but not quite (and also for the next frames there is
> no closeness).
>
> Sorry if this became a long mail, but I need to understand this in order to
> be able to progress with my research.
>
> Thanks
>
>
> ------------------------------
>
> Message: 5
> Date: Tue, 13 Nov 2012 16:54:47 -0500
> From: Justin Lemkul <jalemkul at vt.edu>
> Subject: Re: [gmx-users] Re: Diatomic in MeCN NPT (NH and PR)
>         simulation      segfaults after 1 us
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <50A2C1A7.7050407 at vt.edu>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
>
>
> On 11/13/12 2:56 PM, benjfitz wrote:
> > Are there any other tests I should run to diagnose the problem?
> >
>
> I doubt there's anything that can be done that will be particularly
> useful.  You
> could compile in debugging mode and try to do a backtrace when the problem
> occurs, but if you have a means to run stable simulations (500-ns
> intervals), I
> don't know whether that's worth it to you or not.  I doubt any of the
> developers
> are going to want to run microsecond simulations to try to track it down.
>
> -Justin
>
> --
> ========================================
>
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> ========================================
>
>
> ------------------------------
>
> Message: 6
> Date: Tue, 13 Nov 2012 23:27:28 -0800 (PST)
> From: Raj <princecrusial05 at gmail.com>
> Subject: [gmx-users] Re: hydrophobic contacts
> To: gmx-users at gromacs.org
> Message-ID: <1352878048599-5002931.post at n6.nabble.com>
> Content-Type: text/plain; charset=us-ascii
>
> Hi all,
>
> can some one tel me how can i prepare a index file specifying the
> hydrophobic atoms along for measuring the hydrophobic contacts in the
> systems alone.
>
>
>
> --
> View this message in context:
> http://gromacs.5086.n6.nabble.com/hydrophobic-contacts-tp4998153p5002931.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
>
>
> ------------------------------
>
> --
> gmx-users mailing list
> gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>
> End of gmx-users Digest, Vol 103, Issue 61
> ******************************************
>