[gmx-users] pull code with defined negative relative displacements

[chris.neale at utoronto.ca]

Hello, I am re-running some of our gromacs 3 simulations using gromacs
4, and as far as I can tell the gromacs 4 pull code, while very nicely
enhanced from gromacs 3, has also lost some functionality.

I am calculating the PMF of a peptide across a bilayer and, to  
simplify the issue, what I can't figure out how to do with gromacs 4  
is to pull group1 to a position 1 nm more negative along the z-axis  
than group 0 (the bilayer in this case) for some replicas and to 1 nm  
more positive along the z-axis than group 0 in others. I'll focus on  
the negative displacement here as it is the one that is giving me  
problems.

This used to be possible with the following gromacs 3 pull.ppa options:

runtype   = umbrella
reftype   = com
pulldim   = N N Y
reference_group = group0
group_1   = group1
K1        = 500
Pos1      = 0 0 -1.0

Sure, I could start group1 in the correct negative-z displacement from
group0 and use pull_init1=+1.0, but this will not work when, for
instance, I want to restrain it to -0.1 nm (using pull_init1=+0.1),
where the sampling will infrequently jump back and forth about z=0.

Just to be sure, I tried for gromacs 4 are the following pull code  
.mdp options:

pull                     = umbrella
pull_geometry            = distance
pull_dim                 = N N Y
pull_start               = no
pull_ngroups             = 1
pull_group0              = group0
pull_group1              = group1
pull_init1               = -1.0
pull_rate1               = 0
pull_k1                  = 500.0

where mdrun complains:
"Pull reference distance for group 1 is negative (-1.000000)"

and it is pretty obvious why this doesn't work since it is asking for
a negative displacement. Nevertheless, I tried it and pull_init1
appears to get set to zero.

I also attempted the following.

pull                     = umbrella
pull_geometry            = distance
pull_dim                 = N N Y
pull_start               = no
pull_ngroups             = 1
pull_group0              = group0
pull_group1              = group1
pull_init1               = 1.0
pull_rate1               = 0
pull_k1                  = 500.0
pull_vec1                = 0 0 -1

where I would then use
pull_init1 = 1.0
pull_vec1  = 0 0 1
for the positive side of the bilayer.

However, when I look at the forces, I am getting the negative of what  
I should get when pull_vec1 = 0 0 -1.

coord.xvg:
610.0000	4.9343	-1.02019
660.0001	4.91454	-0.949747

force.xvg:
610.0000	-10.0932
660.0001	25.1265

Although I am getting exactly what I should get when pull_vec1 = 0 0 1  
(for intended positive displacements).

coord.xvg:
660.0001	4.9014	1.16304

force.xvg:
660.0001	-81.5201

Any ideas are greatly appreciated. I can probably mod the code for my
needs, but a standard gromacs binary is always preferable.

Thank you,
Chris.