[gmx-users] Re: Umbrella sampling question

Thu Nov 15 19:44:22 CET 2012

Hi Erik,

Thank you for your answer.

I see your point now. Went and had a look in gmx_wham.c to see how things
are calculated, and that makes sense.

I was looking for an easy way of relating different parts of the resulting
PMF to my original starting frames, as a means to understand exactly what
the PMF is telling me, and I thought that by mapping distances in each
frame to entries in the profile.xvg-file, I could accomplish that. I see
that's not the case.

I am still learning Umbrella sampling and WHAM, and some of it I still do
not get fully. For example, how best (or easiest) to relate (a) events and
structures in individual frames (or in the pulling simulation) to (b) the
pull force vs time plot and (c) the PMF-plot.

Thanks
/PK

2012/11/14 Erik Marklund <erikm at xray.bmc.uu.se>

> Hi,
>
> See below.
>
> 14 nov 2012 kl. 15.06 skrev Gmx QA:
>
> > 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.
> >
>
> profile.xvg doesn't contain distances vs time, so 1.761971 has nothing to
> do with the first frame. Nor is it (very) close to the average distance,
> for the file is energy vs. your reaction coordinate. In your simulations
> 1.761971 is the lowest distance ever encountered and g_wham sets that as a
> reference point for the rest of the free energy profile, hence the 0 in the
> next column. If you inspect the stdout from g_wham you may see how this is
> the case.
>
> Best,
>
> Erik
>
> >
> > 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
> >> ******************************************
> >>
> > --
> > 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!
> > * Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-users-request at gromacs.org.
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> -----------------------------------------------
> Erik Marklund, PhD
> Dept. of Cell and Molecular Biology, Uppsala University.
> Husargatan 3, Box 596,    75124 Uppsala, Sweden
> phone:    +46 18 471 6688        fax: +46 18 511 755
> erikm at xray.bmc.uu.se
> http://www2.icm.uu.se/molbio/elflab/index.html
>
> --
> 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!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>