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

[Thomas Schlesier]

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