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 <[email protected]>:
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
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
