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

[Hadas Leonov]

Dear GMX users,

I am trying to calculate the PMF of a drug inside a membrane protein
which forms a channel. The reaction coordinate is the z-axis. 
In order to do so, I would like to use Umbrella sampling of the drug in
many positions along the z-axis, followed by WHAM. That should give me
the PMF for each position along the pore, inside-out, right?

I am having trouble defining the parameters for the simulations. 
As I understand, I should generate many starting structures, in which
the position of the drug differ by, say, 0.5 Angstrom apart on the z
axis. So let's say I have 30 starting structures (protein + drug + lipid
bilayer + water).

Then Umbrella sampling should be performed on each of those structures,
so I should also have 30 files of pull.ppa. The difference between them
is the position to which the pulled group (drug) is restrained to. Does
that sound reasonable so far?

This is how my pull.ppa looks like: 
--
verbose                  = no
runtype                  = umbrella
ngroups                  = 1
group_1                  = AMAN
;reference_group          = 
weights_1                = 
reference_weights        = 
reftype                  = com
reflag                   = 0
pulldim                  = Y Y N
; DYNAMIC REFERENCE GROUP OPTIONS
r                        = 0
rc                       = 0
update                   = 1

; UMBRELLA SAMPLING OPTIONS       
K1                       = 1000
Pos1                     = 33.482  31.225  33.173
--

K1 is the force constant acting on the pulled group in order to restrain
it, right? how do I determine the best value for it?
For Pos1 - is that OK to input the COM (center of mass) of the ligand,
as it is in each starting structure? Since this is the position I would
like to restrain it to. 

I did not input a reference group, since I'm not sure if the reference
should be the protein or if it'll be better to work in absolute
coordinates.

Pulldim = Y Y N, since for each starting structure, the point is to see
if it moves along x,y due to the forces acting on it by the protein. 
As I understand, since the ligand's z position is different in each
starting structure, we can analyze it with WHAM and calculate the PMF.
Does that sound right, or am I completely wrong?

The result of a simulation with the aforementioned parameters is that
the ligand escapes the protein *into* the lipid environment, which
doesn't make sense. When I use K1=10, the drug does not escape the pore.
So I'm confused here. I thought K1=1000 was supposed to put more force
on the drug to stay in its original position. 

Any help and ideas regarding why this is wrong?

Thanks in advance,
Hadas.