Am 18.04.2012 12:00, schrieb gmx-users-requ...@gromacs.org:
Send gmx-users mailing list submissions to
gmx-users@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-requ...@gromacs.org
You can reach the person managing the list at
gmx-users-ow...@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<aka...@sbc.su.se>
Subject: Re: [gmx-users] g_wham problem with negative COM differences
To: gmx-users@gromacs.org
Message-ID:
<calxandpfxzgyhxdvbsqmsrfgy7-man2f_-vdkbpkrgtzq-t...@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"<jalem...@vt.edu<mailto:jalem...@vt.edu
Subject: Re: [gmx-users] g_wham problem with negative COM
differences
To: Discussion list for GROMACS users<gmx-users@gromacs.org
<mailto:gmx-users@gromacs.org>>
Message-ID:<4f857b2d.3050...@vt.edu
<mailto:4f857b2d.3050...@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'.
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
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
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 gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
