[gmx-users] Re: gmx-users Digest, Vol 73, Issue 188

[Aykut Erbas]

Hi,

Thanks Thomas and Justine for your help


Thomas Schlesier wrote:
> >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...
> >
>
> Yes the problem lies in the simulation setup (or so i think, because i 
> don't know how big your simulation box is and how far you want pull 
> the molecule). But if it is a problem because of pbc the following 
> could help.
> If you have two molecules which are 1nm apart from each other and you 
> want to pull one of it 3nm away, your box length in that direction 
> must be longer then 2 * 4nm, because else you run into problems with pbc.
My system size 8x4x7nm^3, and the size of the peptide is 4.5 nm
The problem lies here I need really large displacements far more larger 
than even the box length  (like I did many times on Gromacs 3.3), in 
this case only way is to increase the system size  which does not look 
reasonable...
 I think currently I will stick to G3.3 and forget about  
periodic_molecules option :)
or I can set pull_pbcatom to an atom at the very edge of the surface  
and gain some more displacement

thanks for your help and interest
A

>
> To get a better feeling about this, take a simple test system:
> Two water molecules in vaccum and pull the away from each other. 
> Center them in the box and them the boxlength apart. Do this one time 
> with pbc and one time without pbc. Then look at the difference in the 
> output of pullf and pullx.
>
> Another thing:
> About the negative force:
> If you pull in the direction of a box vector the force will be 
> positive, if you pull along the other direction (so in the 'minus' 
> direction) the force will be negative.
> To be sure you could use the above test system and pull the other 
> water molecule away. Since they are vaccum simulations, they are 
> really fast, so it doesn't hurt to do it :)
>
> Greetings
> Thomas
>
>
>
> 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