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

chris.neale at utoronto.ca chris.neale at utoronto.ca
Tue Aug 31 01:41:58 CEST 2010


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
>
>
>






More information about the gromacs.org_gmx-users mailing list