Long-range anisotropic dispersion correction

[jchodera]

I just realized we could include the anisotropic long-range dispersion correction in a very simple and inexpensive way via HMC: When we run the dynamics part of MCMC for the MCMCSampler, we always evaluate the reduced potential we want to sample from using a very large cutoff (e.g. just slightly smaller than half the smallest box dimension) when using explicit solvent, but always use a much shorter cutoff (e.g. 9*angstroms) when running dynamics. We can then accept/reject with HMC.

[pgrinaway]

but always use a much shorter cutoff (e.g. 9*angstroms) when running dynamics.

You mean in the NCMC steps?

[jchodera]

but always use a much shorter cutoff (e.g. 9*angstroms) when running dynamics.

For MCMCSampler, let's suppose we alternate between

• update the configuration with some Metropolized MD
• update the discrete state index via the RJMC/NCMC stuff

Suppose we want each update (configuration, state index) to leave the sampler in a "true" (x,k) state sampled from the posterior.

For updating the configuration, we could use HMC:

• Assign velocities from the Maxwell-Boltzmann distribution and constrain positions and velocities
• Compute initial_total_energy with cutoff_long, equal to slightly less than half the smallest box width
• Change the cutoff to cutoff_short = 9A tuned for efficiency
• Propagate dynamics with velocity Verlet for nsteps using timestep, where timestep is tuned for high acceptance
• Change the cutoff to cutoff_long and evaluate final_total_energy
• Accept or reject the whole move with delta_energy = final_total_energy - initial_total_energy in the Metropolis criterion

We might instead choose to use a GHMC for the individual timesteps, in which case we would do this:

• Compute initial_potential_long for cutoff_long
• Change the cutoff to cutoff_short = 9A tuned for efficiency
• Change the cutoff to cutoff_short and evaluate initial_potential_short
• Propagate dynamics with GHMC for nsteps using timestep, where timestep is tuned for high GHMC acceptance
• Compute final_potential_short for cutoff_short
• Change the cutoff to cutoff_long and evaluate final_potential_long
• Accept or reject the whole move with delta_energy = (initial_potential_short - initial_potential_long) + (final_potential_long - final_potential_short) in the Metropolis criterion, discarding the whole sequence of GHMC steps if this is not accepted

[pgrinaway]

Ok, thanks!

[jchodera]

This might have interactions with #77

[jchodera]

Maybe we can have a quick huddle on Monday to figure out what the best way to handle this is (or whether we should punt on it for the first paper)?

[jchodera]

There are several options to consider:

• Ignore it for the first paper. This will cause a systematic deviation from experiment that depends on the number of atoms in the ligand.
• Use the enlarged cutoff in the outermost sampling loops, but a reduced cutoff in MCMC, NCMC, and geometry sampling. This will incur some extra bookkeeping overhead, but will allow us to optimize the smaller cutoff for efficiency. There are currently some issues with enlarging reaction field cutoffs with OpenMM's atomic cutoffs, so we'd either need to stick to PME or just enlarge the Lennard-Jones cutoff (which requires even more complex bookkeeping).
• Just use a big cutoff throughout. If sampling is slow compared to the overhead of other things (e.g. geometry), this may not matter much.

[pgrinaway]

I am inclined to prefer those options in the order 3, 2, 1. My suspicion is that the slowdown might be overwhelmed by the slowness of other components, so it's worth at least taking a look at that option, since it's simple.

Barring that, I think we should try 2, because presumably the reference calculations will have this correction, and I don't think we want to present something that is systematically off.

[jchodera]

We spoke about this in person, and decided we might punt on this until a subsequent paper. We can us the same cutoff (e.g. 9A) for both relative FEP and Perses in the first paper.

[pgrinaway]

@jchodera what should we do about this issue, since we are reweighting to the non-alchemical endpoints anyway?

[jchodera]

We just need to enlarge the nonbonded cutoff for the non-alchemical endpoints. We can enlarge to a little less than half the box size.

[jchodera]

This is correctly handled by unsampled_endstates in our sampler (though omitted the Folding@home calculations we use right now).

@dominicrufa @zhang-ivy : Do you know if the nonequilibrium switching also expands the cutoff at the equilibrium endpoints to include in estimating the reweighted free energies?

[zhang-ivy]

Do you know if the nonequilibrium switching also expands the cutoff at the equilibrium endpoints to include in estimating the reweighted free energies?

It doesn't. I'll need @jchodera and @dominicrufa 's input on how to implement this

[jchodera]

Since we haven't been using the dask distributed version much, the best approach for us would be to

1. Extend the multistate samplers in openmmtools to add a NonequilibriumSwitchingSampler, which just simulates equilibrium endstates and collects work values connecting them for use in MBAR. This will be very simple compared to the other samplers.
2. Use the HybridCompatibilityMixin to extend the NonequilibriumSwitchingSampler in a few lines, as with ReplicaExchangeSampler and SAMSSampler.

[ijpulidos]

CC https://github.com/openmm/openmm/issues/3229

[ijpulidos]

CC https://github.com/openmm/openmm/pull/3236