Stable simulations after geometry proposal

[jchodera]

Just wanted to start a thread to discuss the best way to ensure we get stable simulations after the geometry proposal step.

In order for NCMC to work, we'll need to have a stable simulation for the alchemically-modified lambda=0 state after the geometry proposal step for some timestep. If we can't get stable dynamics with lambda=0, NCMC is not going to help us improve acceptance probabilities. In principle, this should be relatively straightforward since we should only be dealing with intramolecular valence terms at this point.

Obviously, the best solution would be to see if we can improve the geometry proposal, but there are always possibilities where we can't control what is going on with more complex intramolecular interactions or steric clashes that will occur when we start switching lambda from 0 to 1.

Ideas @pgrinaway has suggested are:

• Metropolis Monte Carlo. Rejects NaNs (good!), but does not respect constraints and makes unbiased Gaussian proposals in each atomic coordinate (bad).
• GHMC. Rejects NaNs (good!) and uses gradient information to achieve high acceptance probabilities (also good). Timestep could potentially be autotuned, or even more sophisticated schemes could be used (e.g. extra-chance GHMC). We have to change the NCMC engine to accumulate the appropriate form of work for this, however, since it will be the potential energy change accumulated from each step since Metropolization already handles the shadow work.

Do we want to have the flexibility to try multiple integrators inside of NCMCEngine? If so, what should the API for that look like?

I'm happy to take a stab at this once we decide what the appropriate functionality and API would be.

[pgrinaway]

A couple of extra thoughts here:

• I'm not sure if the reason GHMC and Metropolis work here is that they reject NaNs. I think it's not unlikely that neither of them produce NaNs even in their proposals--in the case of GHMC, if the simulation is heading toward NaNland, it might reject before going there. I'm not sure if the reason here really matters, though.
• Couldn't we get the Metropolis integrator to respect contraints? It might not be worth it though, since as you point out the proposals are kind of crappy.
• Would it be reasonable to simply perform the entire simulation, NCMC and propagation of the coordinates given a compound, with GHMC? It seems to be very well-behaved.
• If we decide to go with GHMC, I think it's worth exploring some sophisticated schemes for tuning the GHMC timestep.
• I really like the idea of an NCMCEngine that can take different integration schemes, since I think this will become very useful for tuning, as well as providing an easy way to benchmark algorithms.

[jchodera]

I'm not sure if the reason GHMC and Metropolis work here is that they reject NaNs. I think it's not unlikely that neither of them produce NaNs even in their proposals--in the case of GHMC, if the simulation is heading toward NaNland, it might reject before going there. I'm not sure if the reason here really matters, though.

For GHMC, you could turn off the Monte Carlo step and see if it explodes. If not, it's the thermostat that is helping.

For both GHMC and MMC, we can add a counter that counts how many times NaNs were rejected.

Couldn't we get the Metropolis integrator to respect contraints? It might not be worth it though, since as you point out the proposals are kind of crappy.

This would be super hard because you either need to make proposals in bond-angle-torsion variables and compute Jacobians for the M-H acceptance criteria, or you need to rewrite the OpenMM constraint algorithm to give you the Lagrange multipliers and compute the Jacobians from that.

Would it be reasonable to simply perform the entire simulation, NCMC and propagation of the coordinates given a compound, with GHMC? It seems to be very well-behaved.

I think we can use GHMC in NCMC by modifying the acceptance criteria, but we have to compute the correct criteria for this. We probably want to understand why it works so much better than VV in case we are missing something.

If we decide to go with GHMC, I think it's worth exploring some sophisticated schemes for tuning the GHMC timestep.

We can definitely add that to the roadmap.

I really like the idea of an NCMCEngine that can take different integration schemes, since I think this will become very useful for tuning, as well as providing an easy way to benchmark algorithms.

Good point. What would the API look like? Which schemes should it support?

[pgrinaway]

For GHMC, you could turn off the Monte Carlo step and see if it explodes. If not, it's the thermostat that is helping.

I think the Monte Carlo step is really what's helping here, since it still explodes when I remove that bit.

This would be super hard because you either need to make proposals in bond-angle-torsion variables and compute Jacobians for the M-H acceptance criteria, or you need to rewrite the OpenMM constraint algorithm to give you the Lagrange multipliers and compute the Jacobians from that.

Ok, nevermind then.

I think we can use GHMC in NCMC by modifying the acceptance criteria, but we have to compute the correct criteria for this. We probably want to understand why it works so much better than VV in case we are missing something.

I agree, I'd like to know as well. Besides the presence of an extra energy computation, is there any reason we should not prefer GHMC or its variants for the propagation of the system (even outside of NCMC)?

Good point. What would the API look like? Which schemes should it support?

I will open another issue for brainstorming this.

[jchodera]

GHMC will be much more expensive than VV for propagating NCMC switches. Maybe 2-10 times per timestep depending on the overhead of energy computation, extra force computation, the global kinetic energy sums, and the accept/reject step.

[jchodera]

It could still end up being more efficient to use GHMC, of course, but from a cost per timestep view, VV is most efficient.

[pgrinaway]

I think we're good here. Can re-open later if necessary.