Holly,
> Thanks for your answer. Actually my fluctuations are worse than you mentioned. Here is what I have
> for fluctuations:
>
> At lambda = 0.01592, 4.7916
> At clambda = 0.08198, 4.9031
> At clambda = 0.19331, 5.6905
> At clambda = 0.33787, 7.5437
> At clambda = 0.5, 11.4195
> At clambda = 0.66213, 24.0413
> At clambda = 0.80669, 48.7498
> At clambda = 0.91802, 133.4857
> At clambda = 0.98408, the program won't run - I think because velocities get too large
>
> So I will have to throw out the values from clambda = 0.91802 on.
> I had previously thought that the RMS fluctuations give the error, but it sounds like this isn't right
> from what you wrote below. I'm not sure how you are estimating errors from the RMS fluctuations.
In terms of the error analysis, the RMS fluctuations are just the
standard deviation of dV/dlambda, basically. These are *not* your
uncertainty. To convince yourself of this, just realize that your
uncertainty should go down as the number of samples of dV/dlambda
increases, but since the RMS fluctuations are just a property of the
distribution of dV/dlambda, these should remain unchanged once you
have enough sampling. So, really, your uncertainty at each lambda
value is (RMS)/sqrt(N/g), where N is your number of samples and g is
the "statistical inefficiency" (always 1 or larger) that measures,
essentially, how many correlated snapshots there are before every
uncorrelated snapshot. In other words, you really don't have N
independent samples of dV/dlambda, so it's incorrect to estimate the
uncertainty by dividing by sqrt(N). Rather, you have N samples of
dV/dlambda, but because the system changes configurations only
relatively slowly (usually), consecutive samples are related, on a
timescale equal to the autocorrelation time (which is why Thomas was
talking about the correlation time).
So, the two basic ways to correctly do error analysis are either (a)
compute the autocorrelation time/statistical inefficiency in order to
use an expression like that above, or (b) only sample dV/dlambda at a
time interval you know to be longer than the autocorrelation time, so
that your samples are uncorrelated and you can just use the sqrt(N)
factor. Of course, usually you don't know the autocorrelation time in
advance so you need to do (a).
> I am really only intersted in the free energy change associated with turning off the WCA attractive
> potential, and so maybe I can do OK with just extrapolating the numbers I have for clambda less than
> or equal to 0.80669 down to zero at clambda = 1, considering that up until clambda = 0.80669 the
> repulsive forces are all still likely on.
In my opinion (and that expressed in our preprint and the Beutler and
van Gunsteren paper I mentioned) the need to extrapolate to the
endpoint is one of the arguments in favor of soft core. But I guess
that doesn't help you very much.
Best wishes,
David
> Best regards,
> Holly
>
>
>
> > Hi,
> >
> > sorry, I am only transiently accesing the net right now and it looks like
> > I deleted half of this discussion before I realised that its on my topic.
> >
> > I think David and Dave have pointed out the important points already,
> > menaing that the fluctuations you see in turning of a vdw-particle
> > (disappearing an atom) are absolutely to be expected. Without klambda>1
> > you have no hope of ever getting a free energy curve that is suitable for
> > numerical integration.
> >
> > The recent paper to which David pointed you indeed discusses precisely
> > that and compares several different scaling schemes, as does the van
> > Gunsteren paper mentioned earlier. From my experience on the solvation
> > free energies, using separation shifted 'softcore' potentials looks like
> > the best option you can have in terms of curve shape that is good for
> > integration and dvdl-histograms that make it likely to get a sound average
> > without doing huge amount of sampling. However, at least for the simple
> > small molecule solvation free energies I looked at, the klambda=6 option
> > is not much worse and should easily give you good results if you dont have
> > access to softcore potentials. I am not sure if this is still the case if
> > a more complex conformational space like a ligand-protein complex is
> > studied (as David also pointed out already). We will do more on a test
> > system there in the near future, but for now I would expect that you have
> > a good chance at least trying to study your transformation with just the
> > klambda option.
> >
> > If I remember right (and correct me if I dont) you mentioned dvdl-rms of
> > ca 30kcal/mol. While this is more than I saw in a disappearing toluene, it
> > is by no means prohibitively large. If you get dvdl-correlation times of
> > say 100fs on a 1 ns run, then this would correspond to a standard error in
> > the range of < 1kcal/mol and that only for the worst part of your free
> > energy curve. So chances are good that you will get a decent free energy
> > with a reasonable standard error.
> >
> > Kind Regards,
> >
> > Thomas
> >
> > >> minor. Thomas is travelling right now, but maybe when he reads this, he
> > >> can
> > >> post the RMS fluctuations from this studies.
> >
> > Sorry, I dont have my data here, so I'll have to refer you to the preprint
> > of our publication available on David Mobley's homepage.
> >
> > -----------------------------------------------------------------------
> > The AMBER Mail Reflector
> > To post, send mail to amber.scripps.edu
> > To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
>
>
>
> --
> Department of Physics, University of Alberta
> Edmonton CANADA
>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber.scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
>
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Jul 08 2007 - 06:07:25 PDT
