Question on relaxation vs alchemical steps

[davidlmobley]

Can I get an update on the latest protocol recommendations (at least in the saltswap context) in terms of number of perturbation vs propagation kernels? As I understand it the issue is that perturbation remains slower than propagation, so there might be a sweet spot where there is still a lot of perturbation, but more propagation. We are exploring related issues for BLUES, and don't want to reinvent the wheel.

@sgill2 and I ran across a discussion of this in openmm/openmm#1832 (comment), where there was a graph of some (apparently older) data on how efficiency varied as a number of perturbation vs propagation kernels from @jchodera. Peter Eastman made the observation that it looked like what was important was the product of the two (though John noted that this was not quite the case). John mentioned that you guys would be exploring this more since perturbation is slower than propagation so there might be a sweet spot where one does more propagation than might otherwise be expected.

(cc @Limn1 )

[jchodera]

Those issues aren't really relevant if you are (1) using alchemical systems generated with one of our alchemical factories, and (2) disable the long-range dispersion calculation for the alchemically-modified atoms in the alchemical factory with the disable_alchemical_dispersion_correction=True argument when creating your alchemically-modified system.

The issues @gregoryross noted were mainly due to the fact that he was instead (1) calling updateParametersInContext and (2) causing the long-range dispersion correction to be numerically recomputed on the CPU each time.

[jchodera]

I'd love to hear about how you guys are doing with BLUES. I think a few weeks back we had asked if you had looked at the stddev of the protocol work as a function of the number of protocol steps, since this is the main way to examine what is going on with protocol efficiency, but I've never seen any plots.

[jchodera]

Specifically, quantifying this for hopping a ligand around in a box of solvent would be sufficient to diagnose any issues.

@gregoryross also found that continuing to use a barostat during the NCMC protocol was important in getting good acceptance rates as well. Not sure if this applies to your case too.

[davidlmobley]

Ah, OK, so I misunderstood slightly. I thought that the thread was saying there was still a major speed difference even aside from the alchemical dispersion correction.

I'd love to hear about how you guys are doing with BLUES.

We've been working through figuring out how to get reasonable acceptance in the "new world" where we no longer use the buggy alchemy which resulted in better acceptance but also crashes due to what I think was basically noise in the work distributions. Looks like a lot more relaxation (20K steps or more) which we are still sorting through.

I think a few weeks back we had asked if you had looked at the stddev of the protocol work as a function of the number of protocol steps, since this is the main way to examine what is going on with protocol efficiency, but I've never seen any plots

I think I'd missed this suggestion, but I'm checking with Sam.

@gregoryross also found that continuing to use a barostat during the NCMC protocol was important in getting good acceptance rates as well. Not sure if this applies to your case too.

Why? This is very interesting -- we have not pursued that at all. Maybe we should look at it.

Specifically, quantifying this for hopping a ligand around in a box of solvent would be sufficient to diagnose any issues.

We've done a bit of that specific test now and can at least get acceptance; we also see the work distribution shift towards "less unfavorable" as we do more relaxation, but have not looked specifically at the standard deviation.

[davidlmobley]

One other thing I've been interested in is whether it would make sense to do more propagation near the "physical" states rather than at the alchemical intermediate states. At least for some simple problems (like rotating a ligand in-place in a cavity in a protein) I can convince myself this would make sense -- for example if you relax a great deal when it is non-interacting, then you may have the issue where the protein cavity collapses around the ligand, and then you have to push it back out of the way again as you turn the ligand back on. In that case, I can wave my hands and convince myself that relaxing more closer to the "physical" states and less towards the intermediate state might result in better acceptance. However, I'm not convinced it's worth exploring this much nor that it would really improve things much.

[jchodera]

We've been working through figuring out how to get reasonable acceptance in the "new world" where we no longer use the buggy alchemy which resulted in better acceptance but also crashes due to what I think was basically noise in the work distributions. Looks like a lot more relaxation (20K steps or more) which we are still sorting through.

This sounds much more in line with my experience and expectations.

One other thing I've been interested in is whether it would make sense to do more propagation near the "physical" states rather than at the alchemical intermediate states. At least for some simple problems (like rotating a ligand in-place in a cavity in a protein) I can convince myself this would make sense -- for example if you relax a great deal when it is non-interacting, then you may have the issue where the protein cavity collapses around the ligand, and then you have to push it back out of the way again as you turn the ligand back on. In that case, I can wave my hands and convince myself that relaxing more closer to the "physical" states and less towards the intermediate state might result in better acceptance. However, I'm not convinced it's worth exploring this much nor that it would really improve things much.

There are likely ways to optimize the protocols to increase acceptance rates---potentially by orders of magnitude---but the first step is to get to the point where the work stddevs are less than a few kT. The relevant theory is explained here:

http://threeplusone.com/Crooks2007c.pdf
http://threeplusone.com/Sivak2012b.pdf

Why? This is very interesting -- we have not pursued that at all. Maybe we should look at it.

We don't know for sure, but the use of a barostat did have a measurable effect.

[sgill2]

I'd love to hear about how you guys are doing with BLUES. I think a few weeks back we had asked if you had looked at the stddev of the protocol work as a function of the number of protocol steps, since this is the main way to examine what is going on with protocol efficiency, but I've never seen any plots.

@jchodera Here's some histogram plots of the work distributions from NCMC using various relaxations steps for toluene in water, and toluene bound to T4 lysozyme.
The titles of each histogram give the following information, in order: "ncmc steps", "acceptance", "total samples", "stddev", "mean". The x axis is the work in terms of -kBT (so positive is more favorable for acceptance).
toluene in water without random rotations

toluene in water with random rotations

toluene in lysozyme with random rotations

[jchodera]

@sgill2: These plots look exactly how I would expect! Thanks!

You should be able to take this data and plot the <P_accept> and its 95% confidence intervals as a function of the number of steps (perhaps on a log-log plot) so that you can see how the acceptance rates are changing with the number of steps. Your acceptance rates are still pretty low in lysozyme, but there will be ways we can increase those acceptance rates.

There are a lot of ways we can increase the acceptance probabilities. Here are just a few ideas:

• Optimize the nonequilibrium path. If you've stored the work vs step number, we can likely estimate the thermoynamic length as a function of lambda and from this develop a better protocol. The easiest way to compute the thermodynamic metric tensor is to estimate var(du/dlambda) at different values of lambda from equilibrium simulations at fixed lambda, however, so you might consider this too. We can also optimize the alchemical softcore parameters, for example.
• Use some method to increase timestep without increasing error. For example, it could be that hydrogen mass partitioning (HMR) lets us take larger timesteps. We're exploring this right now, quantifying the error for a T4 lysozyme system in explicit solvent.
• Don't propagate all degrees of freedom during the switching. The more degrees of freedom there are, the larger the work variance. If we only let residues near the insert-delete locations move (by setting the masses of everything else to zero), for example, we may be able to greatly reduce the work variance because the system behaves like a smaller system.
• Come up with more clever proposals. If you're hopping between endpoints where we don't already have to form a cavity, the dissipative work will be smaller and the acceptance rates will go up.
• Try multiple proposals in parallel. When MC acceptance rates are low, multiple attempts can be made in parallel and one move can be accepted using a modified Metropolis criteria, such as in waste-recycling Monte Carlo. This can linearly increase small acceptance rates.

[sgill2]

@jchodera Thanks! Those all sound like interesting suggestions. I'm excited to try some of them out and see how they affect the acceptance.

[sgill2]

I've ran the same toluene in water NCMC simulations (with random rotations) as above, but this time using the MonteCarloBarostat. I observed higher acceptance rates using the barostat, but the work distributions were distributed quite differently than the NCMC simulations without a barostat.

Toluene in water (with random rotations)

Toluene in water with a barostat (and with random rotations

Toluene in lysozyme (with random rotations)

Toluene in lysozyme with barostat (and with random rotations)

The main difference is that the std deviation doesn't seem to decrease with increased relaxation as before. I'm in the process of running the toluene/lysozyme systems with a barostat, and the shorter relaxation runs show similar behavior to the toluene/water system with a barostat.

We could rationalize the acceptance increase of the barostat in the toluene/water case (since volume changes could help the water respond to the insertion/deletion of toluene), but we can't reason why the acceptance for the lysozyme case would be affected. From what John said earlier about the work distributions, we're leaning towards thinking that this isn't intended behavior and might be due to a bug in the barostat or in our code.

Do you have any thoughts on this @jchodera?

[jchodera]

The main difference is that the std deviation doesn't seem to decrease with increased relaxation as before.

This is actually really concerning, and I'm wondering if this indicates some sort of bug in our work-tracking scheme. Can you point to the code you're using for these tests?

@gregoryross : Can you take a look at the above? I think it is essential to compute the acceptance rate vs number of steps (on a log-log plot) with and without the barostat for saltswap given the results above.

[sgill2]

I'll try to get up a minimum working example up–hopefully later today–using only openmmtools, so you folks don't have to decipher our BLUES code to look into this (and to help determine whether or not this is a bug in BLUES).

[sgill2]

Here’s a minimal example script for toluene in water. In the zip, the example.py script runs for a specified number of relaxation steps and number of iterations and outputs the protocol work after each iteration to a file. You can graph the output using the graph.py script.

bar_example.zip

[davidlmobley]

@gregoryross -- just wanted to confirm you saw this? If not we can follow up with you by e-mail or Slack...?

[jchodera]

I spoke to @gregoryross today. He is going to check his code (which uses ExternalPerturbationLangevinIntegrator) to see if it has the same behavior.

I'm working on running the bar_example.zip myself to verify we see the same thing, and looking into whether something might be going on with AlchemicalNonequilibriumLangevinIntegrator that would cause this issue.

Tagging @pgrinaway too since this issue will directly impact his current perses work, though he is swamped with immediate tasks just to get perses up and running.