<br><div class="gmail_quote"><div class="gmail_extra"><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">
>>>>><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<br>
>>>>> 2.3 to<br>
>>>>> -2.5.<br>
>>>>> ><br>
>>>>> > My problem is that g_wham plots negative COM difference as<br>
>>>>> 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<br>
>>>>> 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>
>>><br>
>> It seems that it is safest if I simply rerun my all simulations...<br>
>> Then we<br>
>> will see how it really is.<br>
>><br>
>> I'm not sure if I understand correctly, but doesn't Thomas's experience<br>
>> also suggest that for some strange algoritmic reason force get's wrong<br>
>> sign<br>
>> when distance gets negative with distance geometry? Or am I completely<br>
>> 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<br>
>> distances (I discarded all runs around zero), I get results that are<br>
>> compatible with my sign trick. And It also made sense biochemically. But<br>
>> well, proper rerun does not leave any questions.<br>
>><br>
>> Does following pull code seem proper to you:<br>
>><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>
>> Regards,<br>
>> Anni<br>
>><br>
><br>
><br>
> It's not really 'strange algoritmic reason', more an algorithmic<br>
> problem. One could say it so: For 'pull_geometry=distance' the problem<br>
> is isotropic, you pull along the vector connecting the two groups. If<br>
> the two groups switch places -> the direction of the vector changes -><br>
> origin of potential jumps<br>
> With 'pull_geometry=direction' you pull along a defined vector. If both<br>
> groups switch places -> but still vector stays the same -> origin of<br>
> potential stays where it was<br>
><br>
> So in theory 'pull_geometry=direction' should solve the problem with<br>
> umbrella sampling a really small distances,<br>
> *BUT* as far as i know 'g_wham' doesn't like 'pull_geometry=direction'.<br>
><br>
<br>
g_wham has supported "direction" and "direction_periodic" for some time. The<br>
only caveat is that one cannot use pullx.xvg files for the WHAM analysis.<br>
<br>
-Justin<br>
<br>
> One idea (but absolut no idea if i will work):<br>
> Run the problematic windows/distances with 'pull_geometry=direction' and<br>
> make a *.tpr file with 'pull_geometry=distance' and just feed this to<br>
> g_wham. One thing which is important is the 'pull_dim'. If it is 1D you<br>
> are fine. If it is 3D, you must calculate the absolute of each force and<br>
> give it the 'right' sign.<br>
><br>
> One thing to note:<br>
> I think to problem occurs only for a handfull of umbrella-sampling<br>
> windows around a distance of 0. And even there, for the larger (but<br>
> still relative small) distances it becomes less problematic (since the<br>
> switching of both groups happens less often).<br>
><br></blockquote><div><br><br>To remind you, I originally used distance geometry for my pmf calculation even<br>I needed negative 'distances' that distance geometry cannot handle. We discussed<br>possibilities to get the proper pmf curve without rerunning everything.<br>
<br>Now I tested my pmf curve with different recreated tpr-files and other settings:<br><br>1. I can get reasonable curve by procedure suggested by you: I recreated tpr-files<br>with position geometry and rerun simulations only at three zero crossing windows.<br>
Here pullx files must be used as input. I will use this curve.<br><br>2. If I use pullf files, recreated tpr-files (either position or direction<br>geometry) and if zero is crossed, then the profile curve is inverted at negative<br>
'distances'. My sign trick inverts the curve back, and the curve seems <br>essentially fine, though I cannot be 100 % sure if it handles region around zero<br>exactly correctly, as there is of course small differences to above case, where<br>
three windows were rerun. <br><br>3. If I recreate tpr-files with position and distance geometry I get curves <br>that are in practice identical, but not exactly identical. I don't understand<br>the reason why they are not exactly identical. My pull codes are shown below. <br>
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.)<br><br>I personally like Direction geometry, because it is simper than Position<br>
geometry, though I got feeling from your earlier answers, that Position geometry<br>would somehow be conceptually more correct. <br><br>Position geometry:<br><br>pull = umbrella<br>pull_geometry = position<br>
pull_dim = N N Y<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 0 0<br>pull_rate1 = 0.0<br>pull_k1 = 1000 ; kJ mol^-1 nm^-2<br> <br>Direction geometry:<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>The original pull code with Distance geometry:<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 pulled<br>pull_init1 = 0<br>pull_rate1 = 0.0<br>pull_k1 = 1000 ; kJ mol^-1 nm^-2<br>
<br><br>Regards,<br>Anni<br><br>... and thank's for your time.<br> </div></div>-- <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>
</div><br><br>