# [gmx-users] Re: Negative pulling forces, periodic molecules

Thomas Schlesier schlesi at uni-mainz.de
Fri May 28 17:39:14 CEST 2010

```I think you run into problems with pbc. the force acting on your pulled
group, is something like
f=k(vt-x)
where x is the displacement from the original position of the pulled
group. i think gromacs does not measure that distance directly but it
measure the distance between your reference and the pulled group. since
is also knows the distance between the reference group and the spring
(spring = the position to where you want to pull) it can calculate the
displacement x and so the force.
i think your problem is the measurement of the distance between the
reference and the pulled group. if the distance is less then half the
box vectors you're fine, if it is longer (in absolute coordinates) you
get a big problem because pbc will make nasty things: the distance
between the two groups chances sign and then the displacement will be
much greater.
A small sketch:
r refercence, p pulled, postion the where it is pulled is somewhere to
the right
1234567
r p      -> p-r = 2
later
1234567
r   p     -> looks like p-r = 4
but is actually
3211234567
p  r        -> p-r = -3
so the force is not k(vt-4) but k(vt+3) which is way higher...

Greetings
Thomas

>Hi, Thanks for your answer, my pull_pbcatoms look like pull_pbcatom0
>belongs to the surface pull_pbcatom2 belongs to the peptide I will be
>happy if you could be more clear about In principle pull_pbcatom0 can
>be any atom, it should not change a predefined pulling rate thanks A.

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

```