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

[Justin Lemkul]

On 6/23/15 8:49 AM, Laura Tociu wrote:
> 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

Not exactly.  The pull rate is not the speed of the pulled group.  A harmonic 
potential is connected to the pulled group, and the virtual spring is moved at 
the specified pull rate.  The pulled group has a force applied to it in 
response, which is added to all the other forces in the system.  The virtual 
spring is always extended a constant rate in this scenario (called 
"constant-velocity pulling/steered MD").  Whether or not the pulled particle(s) 
move(s) at the same rate depends on the other forces in the system and the 
stiffness of the spring.

> 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

The force at t=0 should be zero.  What is its magnitude?

> 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

The spring is increasing in length, increasing the force on it.

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

Depends on the forces involved and the interactions that evolve through the 
trajectory.

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

In pullx.xvg are the coordinates of the COM of the reference (in any relevant 
dimensions, here just z because of your pull-dim settings) and the displacement 
along the axis(axes) specified (again, just z here).

You may want to look at my tutorial for a bit more of the nuts and bolts.

http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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