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

[Justin Lemkul]

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

> 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

========================================