[gmx-users] umbrella sampling (PMF) position discrepancy

Mon Sep 10 17:54:46 CEST 2012

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>
>