[gmx-users] Re: Re[7]: Flat energy profile in g_wham

[chris.neale at utoronto.ca]

First, you confuse me when you say "pulling direction was along  
Y-axis" since your .mdp file indicates pull_dim = Y Y Y and so you are  
pulling in *all* dimensions.

Your profile doesn't look like mine at all, but then again I took mine  
from a force profile. Now that I can see your forces and  
displacements, I can't find any match between them. For example:
cneale at cneale-desktop:~/work/August2010/integral/u10$ tail pullf_u10.xvg
11982.0000	255.769
11984.0000	332.586
11986.0000	291.367
11988.0000	264.031
11990.0000	330.69
11992.0000	336.757
11994.0000	327.46
11996.0000	475.078
11998.0000	440.766
12000.0000	342.814
cneale at cneale-desktop:~/work/August2010/integral/u10$ tail pullx_u10.xvg
11982.0000	3.05189	4.39188	-1.23643	2.03259	2.62201	3.35563
11984.0000	3.12086	4.37951	-1.63297	1.98567	2.64212	-3.25943
11986.0000	3.13274	4.37699	-1.57303	1.97354	2.64089	-3.32607
11988.0000	3.13754	4.38542	-1.52275	1.94917	2.63444	-3.38366
11990.0000	3.12712	4.38416	-1.61647	1.93608	2.63691	-3.29597
11992.0000	3.11342	4.37694	-1.59735	1.95931	2.61705	-3.28952
11994.0000	3.0957	4.384	-1.5691	1.99226	2.57634	-3.3149
11996.0000	3.08384	4.38684	5.24141	2.00727	2.57168	-3.09875
11998.0000	2.99618	4.41111	5.25558	2.08776	2.53869	-3.1227
12000.0000	2.79826	4.41958	5.22444	2.29862	2.53084	-3.12472


Can you explain to me how these can match up (with an equation please)?


Quoting alexander yakovenko <yakovenko_a at ukr.net>:

> Yep, it was mistake with a sign, is   
> http://picasaweb.google.ca/117558205101329348732/G_wham#5511351478636719042   
> looks like your one? Btw my temperature was 310K.
> The dollars you would lost (as I used only 1 mdp for all runs via   
> script and thus avoiding humans mistakes like with pictures) is   
> probably because complicated landscape - the dissociation is coupled  
>  to (partial) unfolding and relative sub-domains shift (actually  
> this  experiment is aimed to support this hypothesis). I.e. it does  
> not  look too suspicious for me; something closer to manual was  
> expected  though. Anyway I attach logs from all my simulations and  
> the mdp.
> As for box, it has 10x16x7 edges and pulling direction was along   
> Y-axis to bring molecules only 5.5nm apart (DNA and protein of ~2 nm  
>  radii) so there is definitely sufficient distances.
> Anyway the only bug in g_wham is an explanation I can give to   
> ABSOLUTELY flat line of PMF profile.
> As for colors, can we please stick on windows ids instead? I guess   
> you told about w4, w6 and w8 so I attached their logs from mdrun and  
>  for w10 too.
> P.S. I still not understand if there is a way to make g_wham working  
>  for me (suggesting all other calculations were made fairly)?
> Alex
>
>
>   First, I don't think that you converted and debiased correctly. Once
> you have a histogram of the displacement, then you will debias that
> histogram like this (assuming run at 298K):
>
> FC=1000
> T=298
> cat biased.hist |awk '{print
> $1,-0.00831447*'${T}'*log($2)-0.5*'${FC}'*$1*$1}' > debiased.hist
>
> You attached the force file, not the coordinate file. I tried to
> convert like this:
>
> cat pullf_u10.xvg|grep -v '[@|#]' |awk '{print $1,-($2/1000)}'
>
> and then I made a histogram and then I converted to energy and
> debiased, but I got something entirely different than yours.
>
> I would bet dollars to donoughts that either your force constant was
> not 1000 or your sampling is otherwise messed up. Can you post the
> coord.xvg from that same window?
>
> What are your .mdp options? if your cutoffs are shorter than 3 nm then
> your 2 molecules are not going to see each other (unless you have
> PME). This is one reason you might get a relatively flat PMF.
> Nevertheless, a perfectly flat PMF is certainly an error! (where's the
> noise?)
>
> My main concern is that your sampling is incorrect for some reason. To
> see this, look at your original set of probability histograms at
> http://picasaweb.google.ca/117558205101329348732/G_wham#5511244561053627250
>
> See the cyan peak on the left half? now look below it and you'll see a
> green (and a purple) histogram line that is half the height and twice
> the width, but is in nearly the same spot. This is suggestive of
> massive unconvergence, so massive that you have probably done
> something wrong.
>
> Please go back through your runs and see if you might have set
> something up incorrectly. If you still think that it's correct, then
> please *cut and paste* sections of your setup for that cyan peak and
> for the green histogram below it (e.g. pull code options in .mdp file,
> commands for mdrun, etc.) and send all of that to the mailing list.
>
> Quoting alexander yakovenko <yakovenko_a at ukr.net>:
>
>> 0. Thank you all who spent time answering my post (and especially
>> for explanations how to convert bins into energy)!
>> 1. Yes, force is 1000 kJ/mol/nm^2 in all windows
>> 2. z is in nm as default in wham
>> 3. I made profile for w#10 (just in the case I attach forces file).
>> the picture of window profile is there
>> http://picasaweb.google.ca/117558205101329348732/G_wham#5511315785714749378.
>> Thank you once more!
>> P.S. Looks like (some of) you are not surprised with zeros in my
>> wham output, but cant you please explain me what is wrong there?
>> P.P.S. Resending message again because of server error, sorry for
>> possible duplication.  
>> Alex
>>  
>> Alexander, please post the following:
>>
>> 1. force constants for your plot at
>> http://picasaweb.google.ca/117558205101329348732/G_wham#5511287294315721970
>> NOTE: if they are all 1000 as in your first post, then you have a
>> problem in your sampling.
>>
>> 2. units for your x-axis in the plot mentioned in part 1, above (i.e.
>> angstroms or nanometers?)
>>
>> 3. convert one (only one) of your sampling probability histograms to
>> an energy histogram and debias it (without using wham). To convert to
>> energy, take each histogram bin and apply E=-kB*T*ln(P). Now debias by
>> subtracting E_rest=0.5*Fc*(r-r0). You will have noise at each end, but
>> what about the middle part, is it flat or not? Post a picture of this
>> please.
>>
>> Chris
>>