[gmx-users] Re: Umbrella sampling question
Erik Marklund
erikm at xray.bmc.uu.se
Wed Nov 14 21:35:36 CET 2012
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
More information about the gromacs.org_gmx-users
mailing list