[gmx-users] Umbrella sampling of a ligand inside a pore
chris.neale at utoronto.ca
Wed Jun 20 21:04:26 CEST 2007
> 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
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
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
More information about the gromacs.org_gmx-users