[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