[gmx-users] Re: umbrella sampling for two polymer interaction

[Thomas Schlesier]

Sorry, my writing was not really excat.
If you use a reference group, the force / potential acting on the pulled 
group is always relative to the reference group.
If you use 'pull_geometry = distance' the origin potential is always in 
the distance 'pull_init1' along the vector from the reference group to 
the pulled group. (assuming the pulling velocity is zero, else there 
would be an additional contribution).
Imagine an MD step. Both molecules have moved.
1) Calculate new origin of the umbrella potential
2) Calculate the forces
3) Let only these forces act which are allowed from 'pull_dim = X Y Z'
Some cases:
Y Y Y -> Move the pulled group to the origin of umbrella potential
Y N N -> let the force act only along the X-axis
N N N -> let no force act
4) do next MD step and start again at (1)

So if one component is N the pulled group can move freely along this 
axis (relative to the reference group).

Just imagine the step of the pull-code as a step which happens after a 
MD step and remember that both step are independed from each other (ok, 
both use the coordinates from the previous step, but they are n that 
sense independend from each other as two following MD steps in a normal 
simulation).
This means, that in the MD step, both groups can move freely, as there 
is now umbrella potential. But in the pull-code step, we switch the 
potential on.

 > Sorry confused. In AMBER I remember I did once methane-methane
 > interaction,just distance based umbrella sampling. But there I did
 > not provide any direction. So should it not be N N N in Gromacs if I
 > want to allow them freely move in xyz.

With 'pull_gemoetry' you tell GROMACS how it should setup the pulling 
potential. With 'distance' it don't needs an explizit pull_vector, since 
GROMACS uses the vector connecting both groups.



Am 31.05.2013 21:04, schrieb gmx-users-request at gromacs.org:
> Dear Thomas,
>
> Thanks a lot again for great reply. Please clarify this also,
>
>>> >>If one uses only 'Y N N' B would only move along the x-axis due to the
> pull, but could move freely in the yz-plane
>>> >>You never want to use the pull-code and 'pull_dim = N N N'
>>> >>This would mean that there is no force acting between your two groups.
> Then one could have skipped using the pull-code...
>
> So keeping Y N N will allow free movement in yz plane. So if I want A B to
> move freely in xyz but just keep them separated by some distance with
> spring const (just like two balloons tied to each other flying in air).
> Sorry confused. In AMBER I remember I did once methane-methane interaction,
> just distance based umbrella sampling. But there I did not provide any
> direction. So should it not be N N N in Gromacs if I want to allow them
> freely move in xyz.
>
> thanks,
>
>
>
> On Fri, May 31, 2013 at 8:53 PM, Thomas Schlesier<schlesi at uni-mainz.de>wrote:
>
>> >For comments to your questions see below.
>> >
>> >More general: (somewhat longer than i wanted. Hope you find some answers
>> >here)
>> >
>> >Imagine two interacting particles A and B which are alinged to the x-axis.
>> >We take A as the reference group, B as pulled group and put the origin of
>> >the umbrella potential on top of B (pull_start=yes).
>> >Simulation starts -> A and B moves.
>> >Pull-code step: From A we calculate the new position of the umbrella
>> >potential, this is unequal to B, since B moved and our reference point move.
>> >Now we have a force acting on B, 'pull_dim' controls in which directions
>> >the force acts. With 'Y Y Y' B is pulled towards the origin of the umbrella
>> >potential (and with this to the position it should have relative to A).
>> >If one uses only 'Y N N' B would only move along the x-axis due to the
>> >pull, but could move freely in the yz-plane. In the end one would get a
>> >structure where A and B have the right distance in the x-axis but are miles
>> >away from each other in the yz-plane.
>> >
>> >Now imagine we pull B away from A. Since MD simulations normally separete
>> >the movement of the center of mass of the system, it would look like A and
>> >B would move away from a middle point.
>> >
>> >(Interchanging A and B should give the same results).
>> >
>> >
>> >Think in your case (doing umbrella sampling) the mdp-file you suggested
>> >would be most appropiate (with 'pull_start=yes' and 'pull_ngroups=1'). This
>> >gives you a potential which fixes the distance between the two proteins.
>> >One thing you should be aware is that if you restrain the distance in 3d,
>> >you have to account for entropic effects (see also the GROMACS manual). If
>> >you restrain the system only in one direction, these don't arise. Think
>> >this is the reason why one sees some work with umbrella sampling were the
>> >restraint works only in one direction.
>> >
>> >
>> >Am 31.05.2013 17:20, schriebgmx-users-request at gromacs.org:
>> >
>> >  Dear Thomas,
>>> >>
>>> >>Thanks a lot for your time and nice explanation. I was not able to get
>>> >>specially the pull_start flag but now its quite clear.
>>> >>
>>> >>I feel sorry, that should be pull_dim = N N N in my case. Also I will be
>>> >>much thankful if you please can help me to make understand following:
>>> >>
>> >
>> >STOP!!!
>> >You never want to use the pull-code and 'pull_dim = N N N'
>> >This would mean that there is no force acting between your two groups.
>> >Then one could have skipped using the pull-code...
>> >
>> >
>> >
>>> >>1)
>>> >>
>>>>>> >>> >>If you do a pulling simulation, there can be reason for chosing the
>>>>> >>>>
>>>> >>>groups (protein = reference , ligand = pulled group, since we want to
>>> >>pull
>>> >>it away)
>>> >>
>>> >>This indeed is correct but I am not able to get depth of this. I mean to
>>> >>say lets keep ligand as a reference and protein as pulled group then yes
>>> >>it
>>> >>sounds stupid but I am not able to provide a reason myself why we can not
>>> >>keep ligand as reference and pull protein rather !!
>>> >>
>> >
>> >Think this setup should also work. For some simple systems i imagine it
>> >should give identical results.
>> >For complex system i would also think so. But i can't comment on this with
>> >actual expirience. The dimer systems which i investigated were symmetric...
>> >
>> >
>> >
>> >
>>> >>
>>> >>2)
>>> >>
>>> >>  > >3) And also what should be pull_ngroups because if there is no
>>>>> >>>>
>>>>>>>>> >>>>> >> >reference group then it should be 2
>>>>>> >>>>>
>>>>>>> >>>> >>
>>>>> >>>>
>>>>> >>> >Better use a reference group -> pull_ngroups = 1
>>>> >>>
>>> >>You don't want to pull in absolute coordinates, when your system can
>>> >>rotate..
>>> >>
>>> >>I am not able to understand this part. Can you please provide some example
>>> >>so that it makes easier to understand this
>>> >>
>> >
>> >Imagine only a single protein which you want to unfold. In an equilibrium
>> >simulation the protein can freely rotate in the box. If we use the
>> >N-terminus as the reference group and the C-terminus a the pulled group,
>> >the origin of the umbrella potential will always be updated and will
>> >account for movement of the N-terminus (reference group).
>> >If one would pull in absolute coordinates, one would need to give the
>> >position of the umbrella potential in absolute space. The molecule can
>> >move, but the origin of the potential will always stay fixed at one place.
>> >Think in the end this would be equal to an position restraint of said group.
>> >If one would want to restrain the distance of two groups in such a way,
>> >one would need two umbrella potentials. But since these would be equal to
>> >two position restraints, there would be no coupling between the two. I
>> >mean, if both groups move around but would have the same distance it should
>> >be ok since the distance is fine. But both umbrella potential would pull
>> >each group back to the initial position.
>> >
>> >
>> >
>>> >>
>>> >>Much thanks again,
>>> >>
>>> >>
>>> >>regards,
>>> >>Jiom
>>> >>
>>> >>
>>> >>
>>> >>
>>> >>On Fri, May 31, 2013 at 1:21 PM, Thomas Schlesier<schlesi at uni-mainz.de**
>>>> >> >wrote:
>>> >>
>>> >>  >Look also into the manual. But the tutorial is a nice place to start.
>>>>> >>> >For further comments see below:
>>>>> >>> >
>>>>> >>> >
>>>>> >>> >  Dear Lloyd,
>>>> >>>
>>>>>>> >>>> >>
>>>>>>> >>>> >>I have read that but my system is different
>>>>>>> >>>> >>
>>>>>>> >>>> >>regards,
>>>>>>> >>>> >>
>>>>>>> >>>> >>
>>>>>>> >>>> >>On Thu, May 30, 2013 at 8:28 PM, lloyd riggs<lloyd.riggs at gmx.ch>
>>>>> >>>>wrote:
>>>>>>> >>>> >>
>>>>>>> >>>> >>  >Dear Jiom,
>>>>> >>>>
>>>>>>>>>> >>>>> >>> >
>>>>>>>>>>> >>>>>> >>> >Look at justines tutorial, there's example pull .mdp.
>>>>>>>>>>> >>>>>> >>> >
>>>>>>>>>>> >>>>>> >>> >Stephan Watkins
>>>>>>>>>>> >>>>>> >>> >
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>
>>>>>> >>>>>
>>>>>>> >>>> >>
>>>>> >>>>
>>>>> >>> >  >
>>>> >>>
>>>>>>>> >>>> >>>
>>>>>> >>>>>
>>>>>>>>>>> >>>>>> >>> >I want to do Umbrella sampling between two different polymers (A
>>>>>>> >>>>>>and B)
>>>>>>>>>>> >>>>>> >>> >interacting with each other with starting configuration
>>>>>>> >>>>>>separated by
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>some
>>>>>> >>>>>
>>>>>>>>>>> >>>>>> >>> >distance and I am trying to bring them closer.
>>>>>>>>>>> >>>>>> >>> >
>>>>>>>>>>> >>>>>> >>> >I have some queries regarding pull inputs: (this is for to run a
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>umbrella
>>>>>> >>>>>
>>>>>>>>>>> >>>>>> >>> >sampling at some distance)
>>>>>>>>>>> >>>>>> >>> >
>>>>>>>>>>> >>>>>> >>> >pull = umbrella
>>>>>>>>>>> >>>>>> >>> >pull_geometry = distance
>>>>>>>>>>> >>>>>> >>> >pull_dim = Y Y Y
>>>>>>>>>>> >>>>>> >>> >pull_start = ???
>>>>>>>>>>> >>>>>> >>> >pull_ngroups = 2?
>>>>>>>>>>> >>>>>> >>> >pull_group0 = polymer_B
>>>>>>>>>>> >>>>>> >>> >pull_group1 = polymer_A
>>>>>>>>>>> >>>>>> >>> >pull_init1 = 0
>>>>>>>>>>> >>>>>> >>> >pull_rate1 = 0.0
>>>>>>>>>>> >>>>>> >>> >
>>>>>>>>>>> >>>>>> >>> >
>>>>>>>>>>> >>>>>> >>> >please suggest for following:
>>>>>>>>>>> >>>>>> >>> >
>>>>>>>>>>> >>>>>> >>> >1) pull_dim I have set to Y Y Y: Is this correct I do not want
>>>>>>> >>>>>>to make
>>>>>>>>>>> >>>>>> >>> >it interact with some directional vector
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>
>>>>>> >>>>>
>>>>>>> >>>> >>
>>>>> >>>>
>>>>> >>> >I don't really understand the question. If you use 'pull_dim = Y Y Y'
>>>> >>>the
>>>>> >>> >pulling potential acts in all 3 dimensions, if you use 'pull_dim = Y N
>>>> >>>N'
>>>>> >>> >the pulling potential acts only in the X direction and your two groups
>>>> >>>an
>>>>> >>> >move freely in the YZ-plane.
>>>>> >>> >
>>>>> >>> >
>>>>> >>> >  >
>>>> >>>
>>>>>>>>> >>>> >>> >2) Which should be group0 or group1, in other words should I pull
>>>>>>> >>>>>>both
>>>>>>>>>>> >>>>>> >>> >together or how I should decide which one should be reference and
>>>>>>>>>>> >>>>>> >>> >which to be pulled as both are different polymers?
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>
>>>>>> >>>>>
>>>>>>> >>>> >>Depends on what you want to do. Easiest way would be define one
>>>>> >>>>polymer a
>>>>> >>>>
>>>>> >>> >group0/reference group and the other as group1/pulled group. System
>>>>> >>> >shouldn't care about which polymer is which group.
>>>>> >>> >If you do a pulling simulation, there can be reason for chosing the
>>>> >>>groups
>>>>> >>> >(protein = reference , ligand = pulled group, since we want to pull it
>>>> >>>away)
>>>>> >>> >
>>>>> >>> >
>>>>> >>> >  >
>>>> >>>
>>>>>>>>> >>>> >>> >3) And also what should be pull_ngroups because if there is no
>>>>>>>>>>> >>>>>> >>> >reference group then it should be 2
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>
>>>>>> >>>>>
>>>>>>> >>>> >>Better use a reference group -> pull_ngroups = 1
>>>>> >>>>
>>>>> >>> >You don't want to pull in absolute coordinates, when your system can
>>>>> >>> >rotate...
>>>>> >>> >
>>>>> >>> >
>>>>> >>> >  >
>>>> >>>
>>>>>>>>> >>>> >>> >4) I am not able to understand pull_start option with pull_init1.
>>>>>>> >>>>>>In
>>>>>>>>>>> >>>>>> >>> >this case if it is set to yes and 0.0 respectively then does
>>>>>>> >>>>>>that mean
>>>>>>>>>>> >>>>>> >>> >this combination is equivalent to pull_start = No if I just
>>>>>>> >>>>>>assume
>>>>>>>>>>> >>>>>> >>> >pull_init1 does not have any default value (which is 0.0); not
>>>>>>>>>>> >>>>>> >>> >existing
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>
>>>>>> >>>>>
>>>>>>> >>>> >> From the setup which you have written above:
>>>>> >>>>
>>>>> >>> >polymer_B is the reference group. the origin of the pulling potential is
>>>>> >>> >at 'pull_init1' (a length) along the vector which goes from polymer_B to
>>>>> >>> >polymer_A (sine you use pull_geometry = distance).
>>>>> >>> >If you set pull_init1=0 and pull_start=no, polymer_A will crash into
>>>>> >>> >polymer_B (since the origin of the umbrella potential is directly at the
>>>>> >>> >center of mass of polymer_B).
>>>>> >>> >If you set pull_init1=0 and pull_start=yes, GROMACS adds the distance
>>>>> >>> >between polymer_B and A to pull_init1 (-> so pull_init1 is now greater
>>>> >>>than
>>>>> >>> >0.0). Now the origin of the umbrella potential is at the center of mass
>>>> >>>of
>>>>> >>> >polymer_A. -> A is restrainted to a certian distance of B.
>>>>> >>> >
>>>>> >>> >  >
>>>> >>>
>>>>>>>>> >>>> >>> >5) Also finally where are upper and lower bounds defined. pull_k1 =
>>>>>>>>>>> >>>>>> >>> >1000 is harmonic applied to some equilibrium distance value. How
>>>>>>> >>>>>>this
>>>>>>>>>>> >>>>>> >>> >distance is taken by the programme (or it is just the starting
>>>>>>>>>>> >>>>>> >>> >distance taken between two groups) and what are the ?? values
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>
>>>>>> >>>>>
>>>>>>>>>>> >>>>>> >>> >defined. (say in AMBER I define r1,r2,r3,r4; where r2=r3 which is
>>>>>>>>>>> >>>>>> >>> >assumed equilibrium value and r1 is lower and r4 is upper value
>>>>>>> >>>>>>which
>>>>>>>>>>> >>>>>> >>> >defines shape of potential)
>>>>>>> >>>>>>
>>>>>>>>> >>>>> >>>
>>>>>> >>>>>
>>>>>>> >>>> >>The umbrella potential is a simple harmonic potential (so no fancy
>>>>> >>>>stuff
>>>>> >>>>
>>>>> >>> >as in AMBER) with
>>>>> >>> >V = 1/2 k x^2
>>>>> >>> >where x is the violation of the equilibrium distance.
>>>>> >>> >For your setup
>>>>> >>> >V = 1/2 (pull_init1 - distance(B-A))^2
>>>>> >>> >where distance(B-A) means the distance between both polymers.
>>>>> >>> >
>>>>> >>> >
>>>>> >>> >Greetings
>>>>> >>> >Thomas