[gmx-users] Re: Re: g_wham problem with negative COM differences
Anni Kauko
akauko at sbc.su.se
Wed May 2 12:49:35 CEST 2012
> >>>>>
> >>>>>
> >>>>> 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).
> >
>
To remind you, I originally used distance geometry for my pmf calculation
even
I needed negative 'distances' that distance geometry cannot handle. We
discussed
possibilities to get the proper pmf curve without rerunning everything.
Now I tested my pmf curve with different recreated tpr-files and other
settings:
1. I can get reasonable curve by procedure suggested by you: I recreated
tpr-files
with position geometry and rerun simulations only at three zero crossing
windows.
Here pullx files must be used as input. I will use this curve.
2. If I use pullf files, recreated tpr-files (either position or direction
geometry) and if zero is crossed, then the profile curve is inverted at
negative
'distances'. My sign trick inverts the curve back, and the curve seems
essentially fine, though I cannot be 100 % sure if it handles region around
zero
exactly correctly, as there is of course small differences to above case,
where
three windows were rerun.
3. If I recreate tpr-files with position and distance geometry I get curves
that are in practice identical, but not exactly identical. I don't
understand
the reason why they are not exactly identical. My pull codes are shown
below.
Which one you would recommend if I would do the whole thing from beginning?
(I might need to do full rerun for other reasons, so then I will be free to
choose the best pull code.)
I personally like Direction geometry, because it is simper than Position
geometry, though I got feeling from your earlier answers, that Position
geometry
would somehow be conceptually more correct.
Position geometry:
pull = umbrella
pull_geometry = position
pull_dim = N N Y
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 0 0
pull_rate1 = 0.0
pull_k1 = 1000 ; kJ mol^-1 nm^-2
Direction geometry:
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
The original pull code with Distance geometry:
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
Regards,
Anni
... and thank's for your time.
--
************************************************
Anni Kauko, Ph.D.
Post Doctoral Researcher
Structural Bioinformatics Laboratory
Dept. of Biosciences, Biochemistry
Åbo Akademi University
20520 Turku, Finland
phone: +358 (0)2 215 4006
mobile: +358 (0)50-576 8656
************************************************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20120502/56532da7/attachment.html>
More information about the gromacs.org_gmx-users
mailing list