Inaccurate decorrelation times when lambda_coul and lambda_vdw change at different times.

[mrshirts]

We performed a calculation turning off an ion with separate lambda_coul and lambda_vdw, and where we first decoupled the coulomb and then the LJ. When alchemical-analysis calculated uncorrelated data points, some of the lambda_coul = 0 / lambda_LJ changing states appeared to have very long correlation times. However, it was computing the correlation time using the sum of dH/dlambda_coul and dH/dlambda_vdw at all lambdas.

Because the LJ soft core was coming off, the lambda_coul spiked quite a bit at an infrequent interval, leading to long correlation times. However, this is spurious, since the dH/dlambda_coul is not used in the calculation of the free energy change.

Upon further inspection, the fact that the sum of the two dH/dlambdas is used affects the lambda_coul = 0 / lambda_LJ = 0 end as well since the decorrelation time is dominated by very fast fluctuations in dH/dlambda_coul, whereas dH/dlambda_LJ has slower fluctuations (which are being washed out).

One possible solution is to only sum together the decorrelation times the lambdas that are changing, since delta G = delta lambda \dot gradient H. That seems a better solution (it matches the observable, where ), but I wanted to verify that other developers agree before trying to make the change. Any other ideas?

[davidlmobley]

Your thought sounds right to me. Want to go ahead and do a PR?

Thanks.

On Thu, Jul 30, 2015 at 3:57 PM, Michael Shirts notifications@github.com
wrote:

We performed a calculation turning off an ion with separate lambda_coul
and lambda_vdw, and where we first decoupled the coulomb and then the LJ.
When alchemical-analysis calculated uncorrelated data points, some of the
lambda_coul = 0 / lambda_LJ changing states appeared to have very long
correlation times. However, it was computing the correlation time using the
sum of dH/dlambda_coul and dH/dlambda_vdw at all lambdas.

Because the LJ soft core was coming off, the lambda_coul spiked quite a
bit at an infrequent interval, leading to long correlation times. However,
this is spurious, since the dH/dlambda_coul is not used in the calculation
of the free energy change.

Upon further inspection, the fact that the sum of the two dH/dlambdas is
used affects the lambda_coul = 0 / lambda_LJ = 0 end as well since the
decorrelation time is dominated by very fast fluctuations in
dH/dlambda_coul, whereas dH/dlambda_LJ has slower fluctuations (which are
being washed out).

One possible solution is to only sum together the decorrelation times the
lambdas that are changing, since delta G = delta lambda \dot gradient H.
That seems a better solution (it matches the observable, where ), but I
wanted to verify that other developers agree before trying to make the
change. Any other ideas?

—
Reply to this email directly or view it on GitHub
#27.

David Mobley
Associate Professor
Department of Pharmaceutical Sciences
Department of Chemistry
3134B Natural Sciences I
University of California, Irvine
Irvine, CA 92697
dmobley@uci.edu
work (949) 824-6383
cell (949) 385-2436

[jchodera]

Why are you using dH/dlambda to compute correlation times if you are using MBAR to analyze the data?

In YANK, we use only the reduced potential u(x_t) if we are running independent simulations, or the sum of the replica reduced potentials at their own thermodynamic states if we are doing replica exchange.

See this IPython notebook for more details.

[davidlmobley]

We're actually using a whole grid of techniques to analyze in general,
including TI (for all supported codes), BAR, MBAR, etc. It makes sense to
use the correlation time of a single quantity for all of the analysis (as
opposed to computing correlation times on dH/dlambda for TI, and on u(x_t)
for BAR/MBAR, for example), we think. My recollection is that there are a
couple reasons we did this on dH/dlambda:
a) in some instances the correlation times seemed to be longer (not
surprising - the potential can fluctuate for reasons very unrelated to the
free energies being computed, whereas the dH/dlambda doesn't do this as
much)
b) some of the codes we've started supporting may not provide data for
BAR/MBAR yet, whereas they all provide data for TI, so it made more sense
to think of TI as the starting point

Thanks.

On Sat, 1 Aug 2015 at 13:50 John Chodera notifications@github.com wrote:

Why are you using dH/dlambda to compute correlation times if you are using
MBAR to analyze the data?

In YANK, we use only the reduced potential u(x_t) if we are running
independent simulations, or the sum of the replica reduced potentials at
their own thermodynamic states if we are doing replica exchange.

See this IPython notebook
http://nbviewer.ipython.org/github/choderalab/simulation-health-reports/blob/master/examples/yank/YANK%20analysis%20example.ipynb
for more details.

—
Reply to this email directly or view it on GitHub
#27 (comment)
.

[jchodera]

You can test this hypothesis. I'd suggest you do an experiment where you see how predictive different choices of observable for the statistical inefficiency are in predicting the true random variation from experiment to experiment (or among uncorrelated blocks of data). You could compute a Q-Q plot similar to the one we use to validate pymbar.

[mrshirts]

So, there are three ready possibilities: dU/dl (assuming the corresponding dhdl), Delta U, and U itself. dU and dU/dL are in most cases going to be very close to each other, since Delta U \approx dU/dL * Delta L for small values of Delta L, and most of the information use to calculate a free energy between two states depends on the dU to nearby states.

Which leaves U and dU/dL I agree we need to do more analysis of this in general, though I agree with david in what we've seen -- U is very noisy, and it can decorrelate faster that dU because of sloshing of energy in the system in DoF that are not involved in the orthogonal DoF. Since MBAR actually only uses dU, not U itself (assuming all at the same T), it's more natural to use dU or dU/dl.

These are only arguments based on limited data, though; the hypothesis can be tested, although the answer is going to be pretty system dependent, so test on a single system won't really give clear answers.

[mrshirts]

So, if we call the first part of the TIprelim code before the decorrelation code, than the lchange array can be used to identify which components contribute to the free energy. Pavel, is there an issue with separating this code out and using it before the uncorrelate subroutine?

Note that dhdl's should only really be used to decorrelate if we know they should be used.

[davidlmobley]

@FEPanalysis - comments on this?

@mrshirts - this sounds reasonable to me but I do not know the code as well as Pavel. Hopefully he will weigh in.

[FEPanalysis]

Now that we have the '-n' flag, should the issue be closed?

[davidlmobley]

For completeness: The -n flag allows decorrelation based on dhdl or dE.

That said, I actually think this is a separate issue, right @mrshirts ? Specifically:

One possible solution is to only sum together the decorrelation times the lambdas that are changing, since delta G = delta lambda \dot gradient H

I think the issue (at least as Michael explains it) is that right now, when we use dhdl for decorelation we use the total dhdl, which might include the dhdl for components with lambda values which are not actually changing (i.e., dhdl for Coulomb lambdas in a case where we're not changing the Coulomb lambda). So really we should be using dhdl for the components which are changing.

This issue is NOT resolved by using the -n flag, as the problem is really that decorrelating based on dhdl is not behaving correctly, and the -n flag bypasses this code to decorrelate based on energy differences. So the dhdl decorrelation is still incorrect.

[mrshirts]

That is a correct summary, David.

On Tue, Dec 1, 2015 at 6:33 PM, David L. Mobley notifications@github.com
wrote:

For completeness: The -n flag allows decorrelation based on dhdl or dE.

That said, I actually think this is a separate issue, right @mrshirts
https://github.com/mrshirts ? Specifically:

One possible solution is to only sum together the decorrelation times the
lambdas that are changing, since delta G = delta lambda \dot gradient H

I think the issue (at least as Michael explains it) is that right now,
when we use dhdl for decorelation we use the total dhdl, which might
include the dhdl for components with lambda values which are not actually
changing (i.e., dhdl for Coulomb lambdas in a case where we're not changing
the Coulomb lambda). So really we should be using dhdl for the components
which are changing.

This issue is NOT resolved by using the -n flag, as the problem is really
that decorrelating based on dhdl is not behaving correctly, and the -n flag
bypasses this code to decorrelate based on energy differences. So the dhdl
decorrelation is still incorrect.

—
Reply to this email directly or view it on GitHub
#27 (comment)
.

[davidlmobley]

To be concrete here, @FEPanalysis , since you're getting going again: To fix this you just need to decorrelate (when decorrelating based on dhdl) based on delta lambda \dot dhdl rather than based on dhdl - i.e., decorrelate only based on the dhdl components which are changing with lambda.