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

Wed Jun 24 15:20:29 CEST 2015

Dear Justin,

Thanks for the reply! Yeah, I understand how the pulling works now. The
forces at time t=0 are not zero, though. There are huge forces such as -270
or even -1700 kj/mol/nm acting on my pull group at time t=0.  What do you
believe could be the cause of that? And why is there a /nm in that force? I
mean, isn't the force per nm (force constant) of the spring always 1000
kj/mol/nm, but the actual force adopts various different values as time
goes by?

I ran a short simulation with pull_coord1_rate = 0, and when I did that I
got reasonable forces such as 40-60 kj/mol/nm (??) but still not a zero
force at time t=0. Please let me know if this is normal behavior or not.

Thanks,

Laura

On Tue, Jun 23, 2015 at 10:27 PM, Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> 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
>
> ==================================================
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>