Poor overlap between alchemical and nonalchemical state when ligand is charged

[andrrizzi]

While performing the analysis for the SAMPLing simulations, I noticed that the 5 replicates estimates are all over the place. The variability is caused by the unsampled states that we use to take into account the PME reciprocal space.

The free energy for the 5 replicates of OA-G3 when excluding the reweighting step are all within ~0.4kcal/mol, but after reweighting they span [-69kcal/mol, 29kcal/mol]. In this calculations, I’m not expanding the cutoff. I’m reweighting to the same cutoff only to include the PME reciprocal space so this is entirely caused by the electrostatics.

I’ve run the overlap check in test_alchemy on the alchemical system with both direct-space and Coulomb PME treatment and here are the results

Direct-space
DeltaF =      251.564 +-        0.999 kT (141.95 samples, g = 4.9); du mean 296.454 kT stddev 16.165 kT

Coulomb
DeltaF =     -578.504 +-        1.000 kT (1839.62 samples, g = 1.1); du mean -50.193 kT stddev 166.910 kT

These are from 2000 samples spaced 1ps from each other.

I'll soon implement in openmmtools the possibility to alchemically modify charges through NonbondedForce.updateParametersInContext (see choderalab/openmmtools#319). This will work for annihilated, non-softened electrostatics, but I opened this thread to brainstorm ideas on how to address this in the general case.

[andrrizzi]

The solution implemented in choderalab/openmmtools#320 and #880 worked. There's now an issue open on the OpenMM's repo with a proposed syntax that may solve this problem without using NonbondedForce.updateParametersInContext(): openmm/openmm#1702 .