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

[Justin A. Lemkul]

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

This is documented, although perhaps not very well.  If you read the manual 
description for the pull_pbcatom1 parameter, this effect is mentioned in a very 
implicit way.

-Justin

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

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================