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

[Raphael Alhadeff]

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.


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

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