[gmx-users] PMF trails off to infinity.

[Thomas Schlesier]

think you encounter the problem, that you construct your pmf from a 3d 
simulation and project it onto 1d, but do no correction.

For TI (if you constrain the distance in all three directions) the pmf 
is given by
V_pmf(r) = - \int [ F_c + 2/(beta*r) ] dr
with F_c the constraint force and \beta = 1/(k_B*T).

for umbrella sampling one should need the same factor
-2/beta \int 1/r dr
if one restraints the system in 3 directions.

Since 'g_wham' uses no *.tpr it can not know the value of 'pull_dim' one 
should introduce the factor afterwards.

There could also the problem that the dummy-atom interacts with the 
ligand, or other clashes, but i don't think so, because the force looks 
too fine. If the ligand would crash into something one should see a 
greatly increasing force. like around 1000-2500 in your force-profile, 
when the protein still helds the ligand in the binding pocket.



Additional comment:
'pull_geometry = distance'
is a really bad idea for pmf generation if one has an actual distance of 
0 and the system is not isotropic, since distance only knows distances 
and has no ideas about directrions:

place a reference particle at the origin of a box
put a second particle on top of it
pull along x
pull along -x
in both cases you get the distance r=|x|=|-x|
if the system is isotropic, everything is fine
if system is anisotropic you get problems

on the nice side, this problem should affect simulations where one pulls 
a ligand away from a binding pocket, only for the windows which is 
centered the nearest to 0.
for mapping an ion-channel with a reference group in the middle of the 
channel it's (/ it should be) far more worse, (if the system is not 
isotropic).
Somewhere in the archive is a longer explainition, i discussed the topic 
with 2-3 different people...






Am 03.07.2012 11:34, schrieb gmx-users-request at gromacs.org:
> Basically what I'm doing is pulling a ligand out of a protein towards a
> dummy atom, which has a mass of 1 and no charge. I've attached the a
> portion .mdp files for both the smd portion and the umbrella sampling. I
> know that the ligand gets very close possibly crashing into the dummy
> atom. So from what you're saying, I'm thinking this might be the source
> of the problem. - Laura ###### smd.mdp########## ; Generate velocites is
> on at 300 K. gen_vel = yes gen_temp = 300.0 gen_seed = 123456 ;COM
> PULLING ; Pull type: no, umbrella, constraint or constant_force pull =
> umbrella ; Pull geometry: distance, direction, cylinder or position
> pull_geometry = distance ; Select components for the pull vector.
> default: Y Y Y pull_dim = Y Y Y ; Cylinder radius for dynamic reaction
> force groups (nm) pull_r1 = 1 ; Switch from r1 to r0 in case of dynamic
> reaction force pull_r0 = 1.5 pull_constr_tol = 1e-06 pull_start = yes
> pull_nstxout = 10 pull_nstfout = 1 ; Number of pull groups pull_ngroups
> = 1 ; Group name, weight (default all 1), vector, init, rate (nm/ps),
> kJ/(mol*nm^2) pull_group0 =JJJ pull_weights0 = pull_pbcatom0 = 0
> pull_group1 = CPZ pull_weights1 = pull_pbcatom1 = 0 pull_vec1 =
> pull_init1 = pull_rate1 = -0.0006 pull_k1 = 800 pull_kB1 = #######
> um.mdp ############## ;COM PULLING ; Pull type: no, umbrella, constraint
> or constant_force pull = umbrella ; Pull geometry: distance, direction,
> cylinder or position pull_geometry = distance ; Select components for
> the pull vector. default: Y Y Y pull_dim = Y Y Y ; Cylinder radius for
> dynamic reaction force groups (nm) pull_r1 = 1 ; Switch from r1 to r0 in
> case of dynamic reaction force pull_r0 = 1.5 pull_constr_tol = 1e-06
> pull_start = yes pull_nstxout = 100 pull_nstfout = 100 ; Number of pull
> groups pull_ngroups = 1 ; Group name, weight (default all 1), vector,
> init, rate (nm/ps), kJ/(mol*nm^2) pull_group0 =JJJ pull_weights0 =
> pull_pbcatom0 = 0 pull_group1 = CPZ pull_weights1 = pull_pbcatom1 = 0
> pull_vec1 = pull_init1 = 0 pull_rate1 = 0 pull_k1 = 1000 pull_kB1 = On
> 07/02/2012 04:38 PM, Justin A. Lemkul wrote:
>> >
>> >  On 7/2/12 4:30 PM, Laura Kingsley wrote:
>>> >>  Just for clarification, the PMF is read from right to left and the force profile
>>> >>  is read from left to right.
>>> >>
>> >  The dramatic change in the magnitude and sign of the force, coupled with the
>> >  steady increase in PMF, indicates to me that some elements of your system are
>> >  crashing into one another.  In the absence of an accompanying explanation of
>> >  what you're doing (description of system, procedure with .mdp parameters, etc)
>> >  that's the best I can offer.
>> >
>> >  -Justin
>> >
>>> >>  On 07/02/2012 04:27 PM, Kingsley, Laura J wrote:
>>>> >>>  Here is a link to both the PMF profile and the force profile:
>>>> >>>
>>>> >>>  http://s1064.photobucket.com/albums/u370/laurakingsley/?action=view&current=pull_fig.jpg
>>>> >>>
>>>> >>>
>>>> >>>
>>>> >>>  -----Original Message-----
>>>> >>>  From:gmx-users-bounces at gromacs.org  [mailto:gmx-users-bounces at gromacs.org] On
>>>> >>>  Behalf Of Justin A. Lemkul
>>>> >>>  Sent: Monday, July 02, 2012 3:56 PM
>>>> >>>  To: Discussion list for GROMACS users
>>>> >>>  Subject: Re: [gmx-users] PMF trails off to infinity.
>>>> >>>
>>>> >>>
>>>> >>>
>>>> >>>  On 7/2/12 3:53 PM, Laura Kingsley wrote:
>>>>> >>>>  Basically what I'm doing is pulling a ligand toward a dummy atom. I pull the
>>>>> >>>>  ligand until its COM is very close (a few angstroms) from the dummy atom. Could
>>>>> >>>>  this be causing this behavior? What information would help? I'm having a tough
>>>> >>>  Collisions between any elements of your system would explain the problem, if it
>>>> >>>  looks like what I'm picturing in my head:)
>>>> >>>
>>>>> >>>>  time getting the figures attached but may be I could email them directly to you?
>>>>> >>>>  Let me know.
>>>>> >>>>
>>>> >>>  Bullet point #4 here provides a suggestion:
>>>> >>>
>>>> >>>  http://www.gromacs.org/Support/Mailing_Lists#Mailing_List_Etiquette
>>>> >>>
>>>> >>>  -Justin
>>>> >>>
> -- Laura Kingsley