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 [email protected]:
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¤t=pull_fig.jpg
>>>
>>>
>>>
>>> -----Original Message-----
>>> From:[email protected] [mailto:[email protected]]
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
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
