[gmx-users] g_wham problem with negative COM differences
Justin A. Lemkul
jalemkul at vt.edu
Wed Apr 18 14:15:53 CEST 2012
Anni Kauko wrote:
>
> 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> <mailto: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>
> > > <mailto:gmx-users at gromacs.org <mailto:gmx-users at gromacs.org>>>
> > > Message-ID: <4F857B2D.3050100 at vt.edu
> <mailto:4F857B2D.3050100 at vt.edu>
> <mailto: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> <http://vt.edu> | (540)
> 231-9080 <tel:%28540%29%20231-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?
>
Using "distance" geometry implies that the sign of the reference displacement
should not change; it must be always positive or always negative to remain
useful. The "direction" and "position" methods are more flexible, as they allow
you to specify the vector along which the restraint is applied, not just the
generic (x,y,z) components.
> 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.
>
I would suggest rerunning with more flexible parameters.
> 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
>
The "direction" geometry should work, but "position" is the most flexible, in my
mind, if you have a molecule that samples configurations with restraint distance
both above and below the center of the reference group. With "direction" and a
vector of (0,0,1) you are telling grompp/mdrun that the restraint is applied
only in the +z dimension, but it seems to me that you need both positive and
negative vectors. Perhaps the right combination of .mdp settings will indeed
give you what you want.
-Justin
--
========================================
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