[gmx-users] Umbrella sampling of a ligand inside a pore

[Chris Neale]

> What I did was to position the ligand in the center of the protein pore,
> and consecutively move it by 0.5 steps along Z (both up and down until
> it reaches the bulk on both sides). In each step I re-center the ligands
> x,y inside the pore to avoid clashes (using nearby residues). Then, each
> of the generated structure is energy-minimized. Hope it's careful
> enough.
>   
Sounds fine. Although I am not sure what you mean by "move it" If you 
are translating it and then you are going to solvate each system 
independently later then that all sounds proper.
> I'll probably perform trial and error before I move to more advanced
> methods. Do I have to optimize the force constant for each simulation or
> is it OK to use the same for all?
>   
Either is ok (as long as your WHAM implementation allows different Fc 
values (most do, but be sure that it does). Start with the same for all 
and if you are missing sampling in some region then you can add a new 
umbrella or perhaps adjust a force constant. I find it easier to predict 
what will happen when I add another umbrella.
> That's the part I'm confused about.
> The reaction coordinate is definitely the z axis. That's why I have 30
> different starting structures of the ligand along Z inside the pore.
> But when you configure pulldim=N N Y, you enable it to move only in the
> z-axis. So you won't see movements of the ligand in the xy plane inside
> the pore (surely the harmonic restraint of x and y wasn't the optimal
> x,y). So isn't it better to use the default pulldim=Y Y Y in that case?
>   
This is not true: "You enable it to move only in the z-axis" What you do 
with NNY is *force* it to move in the z-axis but the x and y can still 
move it's just that they find their own equilibrium.

> The confusing issue is that it seems we're already performing the
> pulling by generating overlapping starting structures, and the remaining
> question is what's the free energy when the ligand is around that
> specific Z (for each Z), yet still let it move freely within XY.
>   
Exactly, but to do that you need to have pulldim=NNY. Here's why:

The accuracy and precision of computationally-derived properties are 
dependent not only on complete sampling of the relevant conformational 
space, but also on the presence of sufficient transitions between local 
minima on the energy landscape so as to provide information on the 
density of states. This implies that unrestrained computer simulations 
can adequately measure only those processes that occur at periods much 
less than the time of data acquisition. And yet, if the relevant motions 
of largest period may be determined and exhaustively sampled by means of 
a biasing potential, the time of data acquisition required from each 
individual simulation is now related to the motions of the second 
largest period. In this case you have (correctly) selected the z axis as 
the reaction coordinate along which the relevant ion motions of largest 
period will occur. You now complete a variety of simulations at varying 
positions along  z, allowing x and y to reach their own equilibrium.
 
If you're still confused (and even if you aren't) then try some of the 
online talks by David Mobley 
(http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial) 
or Alan Grossfield (http://dasher.wustl.edu/alan/talks/wham_talk.pdf) or 
there are lots of others.

Hope that helps, I am about all tapped out on general suggestions for 
this topic.

Chris.