[gmx-users] PMF trails off to infinity.

[Justin A. Lemkul]

On 7/3/12 2:50 PM, Thomas Schlesier wrote:
> 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.
>

The input to g_wham (at least in modern Gromacs versions) is a list of .tpr 
files (passed to the -it flag) and a list of pullx or pullf files (with -ix or 
-if).  How is g_wham ignorant of these pull settings?

-Justin

> 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
>

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================