openmmtools: Inconsistent energies with GHMC integrator

In version 0.8.2 of openmmtools, the custom integrator GHMCIntegrator is able to output energies before and after an integration step. These global variables are called potential_old and potential_new, respectively. I’ve been comparing these energies with ones extracted via context.getState(getEnergy=True). For normal ‘equilibrium’ MD, the energies taken from the integrator match the energies from getState() excellently. However, when modifying the non-bonded parameters with setParticleParameters() and updateParametersInContext(), the energies calculated with the methods above disagree. (I’m modifying the non-bonded parameters for non-equilibrium candidate Monte Carlo.) Tagging @jchodera, @pgrinaway and @juliebehr in case this is related to their Perses project.

I’ve created some minimal code to highlight this behavior. I’m modifying the LJ parameters of a single water molecule in-between two single steps of the GHMC integrator.

from simtk import openmm, unit
from simtk.openmm import app
from openmmtools import integrators
from openmmtools.testsystems import WaterBox

# Creating a box of water:
size = 20.0*unit.angstrom     # The length of the edges of the water box.
temperature = 300*unit.kelvin
pressure = 1*unit.atmospheres
wbox = WaterBox(box_edge=size,nonbondedMethod=app.PME,cutoff=9*unit.angstrom,ewaldErrorTolerance=1E-6)
ghmc = integrators.GHMCIntegrator(temperature, 1/unit.picosecond, 0.1*unit.femtoseconds)

context = openmm.Context(wbox.system, ghmc)
context.setPositions(wbox.positions)
context.setVelocitiesToTemperature(temperature)

force = wbox.system.getForce(2)       # Non-bonded force.

#      Getting the energy BEFORE the parameters are perturbed
ghmc.step(1)
state = context.getState(getEnergy=True)
potential_old = state.getPotentialEnergy()/unit.kilojoule_per_mole
print('old',potential_old,ghmc.getGlobalVariableByName('potential_new'))

#      Perturbation. Reducing LJ parameters by 10%
force.setParticleParameters(0,charge=-0.834,sigma=0.3150752406575124*0.9,epsilon=0.635968*0.9)
force.updateParametersInContext(context)

#      Getting the energy AFTER the parameters have been perturbed.
state = context.getState(getEnergy=True)
potential_new = state.getPotentialEnergy()/unit.kilojoule_per_mole
ghmc.step(1)    
print('new',potential_new,ghmc.getGlobalVariableByName('potential_new'))

#      Resetting the parameters
force.setParticleParameters(0,charge=-0.834,sigma=0.3150752406575124,epsilon=0.635968)
force.updateParametersInContext(context)

As you should see, the printed energies tend to disagree. The effect is more pronounced if the above 'propagation, perturbation, propagation` dynamics is iterated over. I’ve been testing this on CPUs and have some early indications this is occurring on GPUs too.

Any pointers as to why this is happening would be greatly appreciated! Thanks.