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

[Aykut Erbas]

Hi,

Thanks Thomas
I believe you are right... 
when I used g_traj, I also got this non-sense displacement...
And if you check the trajectory on vmd, you see a gradually increasing 
pulling velocity
What is bothering me more is that the problem is not the pull output but 
(I guess) the simulation itself...

thanks

A.


Thomas Schlesier wrote:
> 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
> > >
> > >
> > >
>
>
>