[gmx-users] Negative pulling forces, periodic molecules

[chris.neale at utoronto.ca]

Read about the .mdp options pull_pbcatom0 and pull_pbcatom1

-- original message --

Hi everyone,

I have a pulling code running on single machine with Gromacs 4.0.7
The systems composes of a surface, a protein  and water
I pull the protein from the terminal group on the surface laterally.

With Gromacs 3.3.3 (on single processors ) the setup works perfectly and
generates correct results for pulling positions compared to g_traj...
To be able use periodic molecules (for the surface), I switched to
gromacs 4.0.7 and I use the "position" option to reproduce the previous
v3 results.
The problem  is that the pullf.xvg gives negative forces...
Normally, the pulling force should reach a pseudo-constant level  but in
my case it just keep going toward minus infinity (unless I cease the run)
Also when I checked position of the pulled group I see that the tangent
of the displacement-time curve (which is th pulling rate) is also
increasing afte.
This is bit irrational since the pulling rate is set constant at the
first place...
For instance, with a pulling rate of 10nm/ns in 1ns you should have
displacement of roughly 10nm.
Untill 0.4 ns everything looks ok but then (when the peptide moves to
the next periodic box) velocity increases and I have a displacement of
more 150nm within 1ns...
Also, I must say that in the trajectory the increasing velocity of the
pulled can be seen clearly with naked-eye.

Here is my mdp file;

integrator               = md
nsteps                   = 1000000
dt                       = 0.002
nstlist                  = 40
nstxout                  = 5000
nstvout                  = 5000
nstfout                  = 5000
nstxtcout                = 2000
nstlog                   = 5000
nstenergy                = 5000
constraints              = hbonds
ns_type                  = grid
coulombtype              = pme
pme_order                = 4
fourierspacing           = 0.12
rlist                    = 0.8
rvdw                     = 0.8
rcoulomb                 = 0.8
energygrps               = DIAM Protein SO
tcoupl                   = berendsen
tc_grps                  = DIAM Protein SOL
tau_t                    = 0.1 0.1 0.1
ref_t                    = 300 300 300
gen_vel                  = no
Pcoupl                   = berendsen
Pcoupltype               = semiisotropic
compressibility          = 2.5E-8 4.5E-5
tau_p                    = 1.0 1.0
ref_p                    = 1.0 1.0
comm_mode                = angular
comm_grps                = DIAM
periodic_molecules       = yes
pbc                      = xyz
;PULLING
pull                     = umbrella
pull_start               = yes
pull_geometry            = position
pull_nstxout             = 100
pull_nstfout             = 100
pull_ngroups             = 1
pull_group0              = DIAM
pull_group1              = AA1
pull_vec1                = 1.0 0.0 0.0
pull_init1               = 0.0 0.0 0.0
pull_rate1               = 0.01
pull_k1                  = 100