<br><div class="gmail_quote"><div class="HOEnZb"><div class="h5"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0pt 0pt 0pt 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Anni Kauko wrote:<br>
> ><br>
> > Date: Wed, 11 Apr 2012 08:38:05 -0400<br>
> > From: "Justin A. Lemkul" <<a href="mailto:jalemkul@vt.edu" target="_blank">jalemkul@vt.edu</a> <mailto:<a href="mailto:jalemkul@vt.edu" target="_blank">jalemkul@vt.edu</a>>><br>
> > Subject: Re: [gmx-users] g_wham problem with negative COM<br>
differences<br>
> > To: Discussion list for GROMACS users <<a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br>
> > <mailto:<a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a>>><br>
> > Message-ID: <<a href="mailto:4F857B2D.3050100@vt.edu" target="_blank">4F857B2D.3050100@vt.edu</a><br>
<mailto:<a href="mailto:4F857B2D.3050100@vt.edu" target="_blank">4F857B2D.3050100@vt.edu</a>>><br>
> > Content-Type: text/plain; charset=ISO-8859-1; format=flowed<br>
> ><br>
> ><br>
> ><br>
> > Anni Kauko wrote:<br>
> > > Hi!<br>
> > ><br>
> > > I try to perform pmf calculations for case where a peptide<br>
shifts<br>
> > > through the membrane. My COM differences should vary from 2.3 to<br>
> > -2.5.<br>
> > ><br>
> > > My problem is that g_wham plots negative COM difference as they<br>
> > would be<br>
> > > positive. In pullx-files the COM differences are treated<br>
correctly<br>
> > > (look below). My peptide is not symmetric, so profile curves<br>
are not<br>
> > > symmetric, so loosing the sign for COM difference screws my<br>
profile<br>
> > > curve completely.<br>
> > ><br>
> > > I did not manage to find any pre-existing answers to this<br>
problem<br>
> > from<br>
> > > internet.<br>
> > ><br>
> > > First datalines from pullx files:<br>
> > > (sorry for strange file names...)<br>
> > ><br>
> > > pull_umbr_0.xvg:<br>
> > > 0.0000 6.26031 2.27369<br>
> > ><br>
> > > pullz_umbr_23.xvg:<br>
> > > 0.0000 6.09702 0.0233141<br>
> > ><br>
> > > pullz_umbr_50.xvg:<br>
> > > 0.0000 6.02097 -2.50088<br>
> > ><br>
> > > g_wham command:<br>
> > > g_wham -b 5000 -it tpr_files.dat -ix pullz_files.dat -o<br>
> > > profile_test.xvg -hist histo_test.xvg -unit kCal<br>
> > ><br>
> > > My pull code:<br>
> > ><br>
> > > pull = umbrella<br>
> > > pull_geometry = distance<br>
> > > pull_dim = N N Y<br>
> > > pull_start = yes<br>
> > > pull_ngroups = 1<br>
> > > pull_group0 = POPC_POPS ; reference group is bilayer<br>
> > > pull_group1 = C-alpha_&_r_92-94 ; group that is actually<br>
pulled<br>
> > > pull_init1 = 0<br>
> > > pull_rate1 = 0.0<br>
> > > pull_k1 = 1000 ; kJ mol-1 nm-2<br>
> > ><br>
> ><br>
> > Your problem stems from the use of "distance" geometry. This<br>
method<br>
> > assumes the<br>
> > sign along the reaction coordinate does not change, i.e. always<br>
> > positive or<br>
> > always negative. If the sign changes, this simple method fails.<br>
> > You should be<br>
> > using something like "position" to allow for a vector to be<br>
> > specified. Perhaps<br>
> > you can reconstruct the PMF by separately analyzing the positive<br>
> > restraint<br>
> > distances and negative restraint distances (note here that<br>
> > "distance" really<br>
> > refers to a vector quantity, and thus it can have a sign), or<br>
> > otherwise create<br>
> > new .tpr files using "position" geometry, though I don't know if<br>
> > g_wham will<br>
> > accept them or not.<br>
> ><br>
> > -Justin<br>
> ><br>
> > --<br>
> > ========================================<br>
> ><br>
> > Justin A. Lemkul<br>
> > Ph.D. Candidate<br>
> > ICTAS Doctoral Scholar<br>
> > MILES-IGERT Trainee<br>
> > Department of Biochemistry<br>
> > Virginia Tech<br>
> > Blacksburg, VA<br>
> > jalemkul[at]<a href="http://vt.edu" target="_blank">vt.edu</a> <<a href="http://vt.edu" target="_blank">http://vt.edu</a>> | <a href="tel:%28540%29%20231-9080" value="+15402319080" target="_blank">(540) 231-9080</a><br>
> > <tel:%28540%29%20231-9080><br>
> > <a href="http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin" target="_blank">http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin</a><br>
> ><br>
> > ========================================<br>
> ><br>
> ><br>
> > Thank's!<br>
> ><br>
> > I managed to solve my g_wham problem by doing two things:<br>
> ><br>
> > 1. New tpr-files with proper pull code for g_wham.<br>
> > 2. I also needed to modify signs of pullf values: If value for pullx<br>
> > distance was negative, I reversed the sign of corresponding pullf<br>
value.<br>
> > I did that by my own script.<br>
> ><br>
> > The new pull code:<br>
> ><br>
> > ; Pull code<br>
> ><br>
> > pull = umbrella<br>
> > pull_geometry = direction<br>
> > pull_vec1 = 0 0 1<br>
> > pull_start = yes<br>
> > pull_ngroups = 1<br>
> > pull_group0 = POPC_POPS ; reference group is bilayer<br>
> > pull_group1 = C-alpha_&_r_92-94 ; group that is actually pulled<br>
> > pull_init1 = 0<br>
> > pull_rate1 = 0.0<br>
> > pull_k1 = 1000 ; kJ mol-1 nm-2<br>
> ><br>
> > I am little bit confuced, why I needed to tweak signes of pullf<br>
values.<br>
> > But like that I got the curve that resembles two half curve made for<br>
> > positive and negative pullx distances separately. That curve also<br>
makes<br>
> > sense from biochemical point of view.<br>
> ><br>
>Such changes do not seem appropriate to me. If you change the sign of<br>
>the<br>
>pulling force, you change the implication of what that value means.<br>
>What<br>
>happens if you run your simulations with the new (more appropriate)<br>
>.mdp file?<br>
>Do the forces have the same magnitude, but opposite sign?<br>
<br>
I don't think that the problem is so easy fixed. I had with another<br>
person about a month ago a lengthy discussion on the list.<br>
<br>
The problem is the following:<br>
If you use 'distance' as pull_geometry the position of the minimum of<br>
the umbrella potential is determined by the vector connecting the ref-<br>
and pulled group, but there is no information about the direction.<br>
<br>
Assume this (1D):<br>
reference particle at origin (0)<br>
minimum of umbrella potential (1) - via pull_init1 = 1<br>
pulled particle at (0.5)<br>
-> force acting of pulled particle 0.5*k<br>
-> everything ok<br>
--------------------------------------------------------<br>
we look later and pulled particle moved to (-0.5)<br>
since we are in the same umbrella sampling windows the umbrella<br>
potential minimum should still be at (1) -> force acting would be 1.5*k<br>
(this would be right if we had direction as the pull_geometry)<br>
BUT:<br>
we have pull_geometry=distance<br>
-> minimum of umbrella potential is now at (-1)<br>
(why is this: from ref seen pull-group is now in the negative direction<br>
-> pull_init tells use that we should go along this vector the distance<br>
1 -> new position is now (-1))<br>
-> force acting is now -0.5*k<br>
which is wrong (meaning physical not sensible)<br>
ok, the force is right (from the algorithm standpoint), but the fact<br>
that our minimum of the umbrella potential jumps around during the<br>
sampling screws everything<br>
--------------------------------------------------------------<br>
<br>
hope you get the idea about the origin of problem<br>
<br>
greetings<br>
thomas<br>
<br clear="all"></blockquote></div><br></div></div>It seems that it is safest if I simply rerun my all simulations... Then we will see how it really is. <br><br>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? <br>
<br>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.<br>
<br>Does following pull code seem proper to you:<div class="im"><br><br> pull = umbrella<br> pull_geometry = direction<br> pull_vec1 = 0 0 1<br> pull_start = yes<br> pull_ngroups = 1<br> pull_group0 = POPC_POPS ; reference group is bilayer<br>
pull_group1 = C-alpha_&_r_92-94 ; group that is actually pulled<br> pull_init1 = 0<br> pull_rate1 = 0.0<br> pull_k1 = 1000 ; kJ mol-1 nm-2<br>
<br><br></div>Regards,<br>Anni<br><br>-- <br>************************************************<br>Anni Kauko, Ph.D.<br>Post Doctoral Researcher<br><br>Structural Bioinformatics Laboratory<br>Dept. of Biosciences, Biochemistry<br>
Åbo Akademi University<br>20520 Turku, Finland<br>phone: <a href="tel:%2B358%20%280%292%20215%204006" value="+35822154006" target="_blank">+358 (0)2 215 4006</a><br>mobile: <a href="tel:%2B358%20%280%2950-576%208656" value="+358505768656" target="_blank">+358 (0)50-576 8656</a><br>
<br>************************************************<br>
</div><br>