[gmx-users] Re: g_wham problem with negative COM differences
Justin A. Lemkul
jalemkul at vt.edu
Wed Apr 18 14:43:51 CEST 2012
Thomas Schlesier wrote:
> Am 18.04.2012 12:00, schrieb gmx-users-request at gromacs.org:
>> Send gmx-users mailing list submissions to
>> gmx-users at gromacs.org
>>
>> To subscribe or unsubscribe via the World Wide Web, visit
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> or, via email, send a message with subject or body 'help' to
>> gmx-users-request at gromacs.org
>>
>> You can reach the person managing the list at
>> gmx-users-owner at gromacs.org
>>
>> When replying, please edit your Subject line so it is more specific
>> than "Re: Contents of gmx-users digest..."
>>
>>
>> Today's Topics:
>>
>> 1. Re: g_wham problem with negative COM differences (Anni Kauko)
>> 2. Error- Atom not found in residue seq nr while adding atom
>> (aiswarya pawar)
>>
>>
>> ----------------------------------------------------------------------
>>
>> Message: 1
>> Date: Wed, 18 Apr 2012 12:19:08 +0300
>> From: Anni Kauko<akauko at sbc.su.se>
>> Subject: Re: [gmx-users] g_wham problem with negative COM differences
>> To: gmx-users at gromacs.org
>> Message-ID:
>> <CALxANdPFxZgYhxdVbsqMSrfgY7-MaN2F_-VDkbpkrgtZQ-tYUw at mail.gmail.com>
>> Content-Type: text/plain; charset="iso-8859-1"
>>
>>> Anni Kauko wrote:
>>>>>
>>>>> Date: Wed, 11 Apr 2012 08:38:05 -0400
>>>>> From: "Justin A. Lemkul"<jalemkul at vt.edu<mailto:jalemkul at vt.edu
>>>>>
>>>>> Subject: Re: [gmx-users] g_wham problem with negative COM
>>> differences
>>>>> To: Discussion list for GROMACS users<gmx-users at gromacs.org
>>>>> <mailto:gmx-users at gromacs.org>>
>>>>> Message-ID:<4F857B2D.3050100 at vt.edu
>>> <mailto:4F857B2D.3050100 at vt.edu>>
>>>>> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>>>>
>>>>>
>>>>>
>>>>> Anni Kauko wrote:
>>>>> > Hi!
>>>>> >
>>>>> > I try to perform pmf calculations for case where a peptide
>>> shifts
>>>>> > through the membrane. My COM differences should vary from
>>>>> 2.3 to
>>>>> -2.5.
>>>>> >
>>>>> > My problem is that g_wham plots negative COM difference as
>>>>> they
>>>>> would be
>>>>> > positive. In pullx-files the COM differences are treated
>>> correctly
>>>>> > (look below). My peptide is not symmetric, so profile curves
>>> are not
>>>>> > symmetric, so loosing the sign for COM difference screws my
>>> profile
>>>>> > curve completely.
>>>>> >
>>>>> > I did not manage to find any pre-existing answers to this
>>> problem
>>>>> from
>>>>> > internet.
>>>>> >
>>>>> > First datalines from pullx files:
>>>>> > (sorry for strange file names...)
>>>>> >
>>>>> > pull_umbr_0.xvg:
>>>>> > 0.0000 6.26031 2.27369
>>>>> >
>>>>> > pullz_umbr_23.xvg:
>>>>> > 0.0000 6.09702 0.0233141
>>>>> >
>>>>> > pullz_umbr_50.xvg:
>>>>> > 0.0000 6.02097 -2.50088
>>>>> >
>>>>> > g_wham command:
>>>>> > g_wham -b 5000 -it tpr_files.dat -ix pullz_files.dat -o
>>>>> > profile_test.xvg -hist histo_test.xvg -unit kCal
>>>>> >
>>>>> > My pull code:
>>>>> >
>>>>> > pull = umbrella
>>>>> > pull_geometry = distance
>>>>> > pull_dim = N N Y
>>>>> > pull_start = yes
>>>>> > pull_ngroups = 1
>>>>> > pull_group0 = POPC_POPS ; reference group is
>>>>> bilayer
>>>>> > pull_group1 = C-alpha_&_r_92-94 ; group that is actually
>>> pulled
>>>>> > pull_init1 = 0
>>>>> > pull_rate1 = 0.0
>>>>> > pull_k1 = 1000 ; kJ mol-1 nm-2
>>>>> >
>>>>>
>>>>> Your problem stems from the use of "distance" geometry. This
>>> method
>>>>> assumes the
>>>>> sign along the reaction coordinate does not change, i.e. always
>>>>> positive or
>>>>> always negative. If the sign changes, this simple method fails.
>>>>> You should be
>>>>> using something like "position" to allow for a vector to be
>>>>> specified. Perhaps
>>>>> you can reconstruct the PMF by separately analyzing the positive
>>>>> restraint
>>>>> distances and negative restraint distances (note here that
>>>>> "distance" really
>>>>> refers to a vector quantity, and thus it can have a sign), or
>>>>> otherwise create
>>>>> new .tpr files using "position" geometry, though I don't know if
>>>>> g_wham will
>>>>> accept them or not.
>>>>>
>>>>> -Justin
>>>>>
>>>>> --
>>>>> ========================================
>>>>>
>>>>> Justin A. Lemkul
>>>>> Ph.D. Candidate
>>>>> ICTAS Doctoral Scholar
>>>>> MILES-IGERT Trainee
>>>>> Department of Biochemistry
>>>>> Virginia Tech
>>>>> Blacksburg, VA
>>>>> jalemkul[at]vt.edu<http://vt.edu> | (540) 231-9080
>>>>> <tel:%28540%29%20231-9080>
>>>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>>>>
>>>>> ========================================
>>>>>
>>>>>
>>>>> Thank's!
>>>>>
>>>>> I managed to solve my g_wham problem by doing two things:
>>>>>
>>>>> 1. New tpr-files with proper pull code for g_wham.
>>>>> 2. I also needed to modify signs of pullf values: If value for pullx
>>>>> distance was negative, I reversed the sign of corresponding pullf
>>> value.
>>>>> I did that by my own script.
>>>>>
>>>>> The new pull code:
>>>>>
>>>>> ; Pull code
>>>>>
>>>>> pull = umbrella
>>>>> pull_geometry = direction
>>>>> pull_vec1 = 0 0 1
>>>>> pull_start = yes
>>>>> pull_ngroups = 1
>>>>> pull_group0 = POPC_POPS ; reference group is bilayer
>>>>> pull_group1 = C-alpha_&_r_92-94 ; group that is actually pulled
>>>>> pull_init1 = 0
>>>>> pull_rate1 = 0.0
>>>>> pull_k1 = 1000 ; kJ mol-1 nm-2
>>>>>
>>>>> I am little bit confuced, why I needed to tweak signes of pullf
>>> values.
>>>>> But like that I got the curve that resembles two half curve made for
>>>>> positive and negative pullx distances separately. That curve also
>>> makes
>>>>> sense from biochemical point of view.
>>>>>
>>>> Such changes do not seem appropriate to me. If you change the sign of
>>>> the
>>>> pulling force, you change the implication of what that value means.
>>>> What
>>>> happens if you run your simulations with the new (more appropriate)
>>>> .mdp file?
>>>> Do the forces have the same magnitude, but opposite sign?
>>>
>>> I don't think that the problem is so easy fixed. I had with another
>>> person about a month ago a lengthy discussion on the list.
>>>
>>> The problem is the following:
>>> If you use 'distance' as pull_geometry the position of the minimum of
>>> the umbrella potential is determined by the vector connecting the ref-
>>> and pulled group, but there is no information about the direction.
>>>
>>> Assume this (1D):
>>> reference particle at origin (0)
>>> minimum of umbrella potential (1) - via pull_init1 = 1
>>> pulled particle at (0.5)
>>> -> force acting of pulled particle 0.5*k
>>> -> everything ok
>>> --------------------------------------------------------
>>> we look later and pulled particle moved to (-0.5)
>>> since we are in the same umbrella sampling windows the umbrella
>>> potential minimum should still be at (1) -> force acting would be 1.5*k
>>> (this would be right if we had direction as the pull_geometry)
>>> BUT:
>>> we have pull_geometry=distance
>>> -> minimum of umbrella potential is now at (-1)
>>> (why is this: from ref seen pull-group is now in the negative direction
>>> -> pull_init tells use that we should go along this vector the distance
>>> 1 -> new position is now (-1))
>>> -> force acting is now -0.5*k
>>> which is wrong (meaning physical not sensible)
>>> ok, the force is right (from the algorithm standpoint), but the fact
>>> that our minimum of the umbrella potential jumps around during the
>>> sampling screws everything
>>> --------------------------------------------------------------
>>>
>>> hope you get the idea about the origin of problem
>>>
>>> greetings
>>> thomas
>>>
>>>
>> It seems that it is safest if I simply rerun my all simulations...
>> Then we
>> will see how it really is.
>>
>> I'm not sure if I understand correctly, but doesn't Thomas's experience
>> also suggest that for some strange algoritmic reason force get's wrong
>> sign
>> when distance gets negative with distance geometry? Or am I completely
>> confuced? If it is so, wouldn't my trick then be justified?
>>
>> At least if I plot two half curves separately for positive and negative
>> distances (I discarded all runs around zero), I get results that are
>> compatible with my sign trick. And It also made sense biochemically. But
>> well, proper rerun does not leave any questions.
>>
>> Does following pull code seem proper to you:
>>
>>
>> pull = umbrella
>> pull_geometry = direction
>> pull_vec1 = 0 0 1
>> pull_start = yes
>> pull_ngroups = 1
>> pull_group0 = POPC_POPS ; reference group is bilayer
>> pull_group1 = C-alpha_&_r_92-94 ; group that is actually pulled
>> pull_init1 = 0
>> pull_rate1 = 0.0
>> pull_k1 = 1000 ; kJ mol-1 nm-2
>>
>>
>> Regards,
>> Anni
>>
>
>
> It's not really 'strange algoritmic reason', more an algorithmic
> problem. One could say it so: For 'pull_geometry=distance' the problem
> is isotropic, you pull along the vector connecting the two groups. If
> the two groups switch places -> the direction of the vector changes ->
> origin of potential jumps
> With 'pull_geometry=direction' you pull along a defined vector. If both
> groups switch places -> but still vector stays the same -> origin of
> potential stays where it was
>
> So in theory 'pull_geometry=direction' should solve the problem with
> umbrella sampling a really small distances,
> *BUT* as far as i know 'g_wham' doesn't like 'pull_geometry=direction'.
>
g_wham has supported "direction" and "direction_periodic" for some time. The
only caveat is that one cannot use pullx.xvg files for the WHAM analysis.
-Justin
> One idea (but absolut no idea if i will work):
> Run the problematic windows/distances with 'pull_geometry=direction' and
> make a *.tpr file with 'pull_geometry=distance' and just feed this to
> g_wham. One thing which is important is the 'pull_dim'. If it is 1D you
> are fine. If it is 3D, you must calculate the absolute of each force and
> give it the 'right' sign.
>
> One thing to note:
> I think to problem occurs only for a handfull of umbrella-sampling
> windows around a distance of 0. And even there, for the larger (but
> still relative small) distances it becomes less problematic (since the
> switching of both groups happens less often).
>
> Greetings
> Thomas
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list