[gmx-users] umbrella sampling (PMF) position discrepancy
Jianguo Li
ljggmx at yahoo.com.sg
Tue Sep 11 03:30:21 CEST 2012
I guess is that g_wham takes the distance from tpr file which calculates the distance using grommp.
It seems that the distances calculated using g_dist and grompp are different, as discussed in this forum about 10 days ago.
-Jianguo
----- Original Message -----
From: Raphael Alhadeff <raphael.alhadeff at mail.huji.ac.il>
To: Discussion list for GROMACS users <gmx-users at gromacs.org>
Cc:
Sent: Monday, 10 September 2012, 23:54
Subject: Re: [gmx-users] umbrella sampling (PMF) position discrepancy
Hi Justin, thank you for you quick reply.
On Mon, Sep 10, 2012 at 6:41 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> On 9/10/12 11:22 AM, Raphael Alhadeff wrote:
>
>> Dear Gromacs users,
>>
>> I've been trying to run an umbrella sampling (for the purpose of pmf)
>> using
>> Gromacs 4.5.5.
>> My system consists of a membrane protein transporter and a Na ion passing
>> through it, from the water bulk on one side of the membrane to the water
>> bulk on the other side. I've used the protein as the reference group and
>> the Na ion as the pull group 1 (I've attached the mdp file below). I have
>> used position geometry because to my understanding it is good for the case
>> where the pulled group crosses the reference's COM. I made 31 frames where
>> the ion is <0.2 nm apart, from -2 nm to +2 (in respect to the COM of the
>> protein). I confirmed these distances using g_traj and g_dist on the
>> starting gro files.
>>
>> As I run the g_wham analysis, using -v, I see that the 'position' of each
>> of my frames are highly different than the distances I've measured in the
>> gro files, not only in value but also relative to each other. The same
>> numbers appear in the pull_x or pull_f (after converting appropriately)
>> files.
>> So what I don't understand is how does Gromacs calculate the values for
>> pull_x (and thus for the 'position' for g_wham). I was under the
>> impression
>> it uses COM distance between group0 and group1, but trying to compare
>> using
>> g_traj or g_dist proved me wrong.
>> I have pasted 5 datapoints as an example below, giving the distance
>> measured by g_dist and the distance that pull_x gives:
>>
>> time(ps) g_dist(z) pull_x(1dz)
>> 0 0.894 -1.817
>> 10 0.857 -1.698
>> 20 0.897 -1.866
>> 30 0.890 -1.913
>> 40 0.966 -1.781
>> 50 0.819 -1.76
>>
>>
>> I've read countless of threads before posting this, and could not find any
>> answer, and will be very appreciative for some light into this.
>>
>>
> The result of g_dist is the positive root of the distance equation. It
> also uses x, y, and z components of the distance, while in your case only
> the z component may be relevant. The dist.xvg file(s) will have each
> component listed after the total distance in subsequent columns.
>
> I understand, and that is what I think I posted in the example. I have the
z component of g_dist and following that my umbrella simulation pull_x 1dz
value (which according to the parameters I gave in the mdp, should to the
best of my knowledge, give the same number) yet the numbers are quite
different.
>
>> I should mention that I am not using positional restraints on the protein.
>> I assume that since the protein is inside the membrane and the ion is much
>> smaller than the protein, the PR is not required, and I wanted the protein
>> to be able to adjust slightly to the movement of the ion (this is a
>> transporter, not a channel).
>> If this is the reason that is causing me trouble I will be happy to have a
>> short explanation on why this makes the distances seem somewhat random.
>>
>> Lastly, I will use this opportunity on the forum to have 2 technical
>> clarifications -
>> -Using pull_init = 0 (or any other number) does not overrun pull_start,
>> rather it adds up, correct?
>>
>
> Correct. The output of grompp prints the distances that will be used, as
> a way to check.
>
>
> -What is the difference between the profile.xvg created (default name) and
>> the pmfintegrated.xvg that is sometimes being created, and how does g_wham
>> decide whether to create one or not?
>>
>>
> Which flag is giving you pmfintegrated.xvg? That's not a default name, so
> without your g_wham command line, we can't guess.
>
> -Justin
>
> That is what I don't understand, I gave no flag for this file, my command
line is simply
g_wham -ix pull_x.dat -it tpr.dat -v
Thank you again,
Raphael
> Thank you very much for the help..
>>
>> Raphael
>>
>> mdp file:
>>
>> title = pmf
>>
>> integrator = md
>> dt = 0.002
>>
>> nsteps = 5000000 ; 10 ns
>> nstcomm = 10
>>
>> ; Output parameters
>> nstxout = 5000 ; every 10 ps
>> nstvout = 5000
>> nstxtcout = 5000 ; every 1 ps
>> nstenergy = 5000
>>
>> ; Bond parameters
>> constraint_algorithm = lincs
>> constraints = all-bonds
>> continuation = yes ; continuing from NPT
>>
>> ; Single-range cutoff scheme
>> nstlist = 5
>> ns_type = grid
>> rlist = 1.2
>> rcoulomb = 1.2
>> rvdw = 1.2
>>
>> ; PME electrostatics parameters
>> coulombtype = PME
>> fourierspacing = 0.16
>> pme_order = 4
>>
>> ; Berendsen temperature coupling is on in two groups
>> Tcoupl = Nose-Hoover
>> tc_grps = Protein POP Sol_Ions
>> tau_t = 0.5 0.5 0.5
>> ref_t = 310 310 310
>>
>> ; Pressure coupling is on
>> Pcoupl = Parrinello-Rahman
>> pcoupltype = semiisotropic
>> tau_p = 2.0 2.0
>> compressibility = 4.5e-5 4.5e-5
>> ref_p = 1.0 1.0
>>
>> ; Generate velocities is off
>> gen_vel = yes
>> gen_seed = -1
>>
>> ; Periodic boundary conditions are on in all directions
>> pbc = xyz
>>
>> ; Long-range dispersion correction
>> DispCorr = EnerPres
>>
>> ; Pull code
>> pull = umbrella
>> pull_geometry = position
>> pull_dim = N N Y
>> pull_start = yes
>> pull_ngroups = 1
>> pull_group0 = Protein
>> pull_group1 = r_4608
>> pull_init1 = 0 0 0
>> pull_vec1 = 0 0 1
>> pull_rate1 = 0.0
>> pull_k1 = 1000
>>
>>
> --
> ==============================**==========
>
> 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<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
>
> ==============================**==========
> --
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
> * Please search the archive at http://www.gromacs.org/**
> Support/Mailing_Lists/Search<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<http://www.gromacs.org/Support/Mailing_Lists>
>
--
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
More information about the gromacs.org_gmx-users
mailing list