[gmx-users] The pullf.xvg and pullx.xvg files

[Laura Tociu]

Dear Gromacs users,

I would like to understand how exactly the pull code works and what exactly
is outputted in the two .xvg files. I am using Gromacs 5.0.2 and I am
pulling an ion through a membrane channel in the xy plane, and the
parameters I am using seem to be working well:

; Pull code
pull    = umbrella
pull-geometry   = cylinder
pull-coord1-vec = 0 0 1
pull-dim        = N N Y
pull-start      = yes
pull-ngroups    = 2
pull-ncoords    = 1
pull-coord1-groups      = 1 2
pull-group1-name        = Protein
pull-group2-name        = CAIon
pull-coord1-rate        = -0.0005
pull-coord1-k           = 1000
pull-r1         = 2
pull-r0         = 3.2
pull-print-reference    = yes

I can see the ion move through the channel by using the above pull code, so
everything should be ok.

However, I am terribly confused by the pullf.xvg and the pullx.xvg files.
To start with what baffles me the most, let me first explain how I
understand the pulling to work. From my understanding at time t = 0, I have
the pull_group, namely my ion, and the reference group, the protein
channel, sitting quite comfortably somewhere in space. Then, at time t = 2
femtoseconds or whatever the time step I am using is, the ion is moved to a
new position based on the regular forces that act on it during MD. However,
this new position will not correspond precisely to the position at which it
should be at 2 fs if it were to travel at the pull rate I indicated in the
.mdp file and in the pull direction I indicated in the .mdp file. So this
is where the umbrella potential comes into play, adding an extra force that
is trying to force the ion in the direction and at the speed I intended it
to go in the .mdp file, and the final location will be one determined by
all forces together. That is to say, the bottom of the "umbrella" potential
is always at the position at which the ion is supposed to be according to
the pull rate, the pull vector and the time t. At first the force due to
the potential will be very small, but over time, if the ion fails to travel
as desired, this force will amount until finally the ion starts to break
through anything that may be blocking it. If this were indeed how it works,
the harmonic force acting on the ion at t=0 should be zero. Nonetheless, in
the pullf.xvg file the force at time t=0 is never zero, but seems to
increase in magnitude the further away the ion is from the COM of the
protein channel (well, some cylindrical part of it anyway). This seems to
indicate that the umbrella potential somehow has its bottom at the COM of
the reference group, which I am not sure I understand the purpose of. For
example in one case where I started with the ion very far away from the
channel, the force at time t=0 was -1700 KJ/mol, and I don't understand why
such a big force would be required to move it at a constant rate towards
the channel, when with the same settings but an ion closer to the channel
the initial force is much smaller (-270 KJ/mol).

I am also unsure about exactly what the "1 cZ" and "1 dZ" columns refer to
in the pullx.xvg file, but I would like to understand the forces first, as
those will get used next in the actual umbrella sampling and the WHAM
analysis.

Any input/clarification would be extremely helpful.

Best,

Laura