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

[Thomas Schlesier]

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, schrieb gmx-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