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

[Thomas Schlesier]

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

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