Parameters reset on reinitialization

[sgill2]

I've been trying to do the following: 1.Take i number of steps with openmmtools. AlchemicalNonequilibriumLangevinIntegrator.

2. After i steps, set the masses of a subset of the system to zero.
3. Call context.reinitialize() so the changes to the system's masses are actually applied to the simulation.
4. Pass in the saved state before reinitialize so after reinitialization the context is the same as before (except the masses).

However when I coded this up I got behavior I didn't expect, namely that after reinitializing the lambda_sterics and lambda_electrostatics parameters were reset to 1.0.

Some example print output related to the parameters

parameters ['lambda_electrostatics', 'lambda_sterics']
parameters before reinitialization [0.0, 0.6666666666666667]
parameters after reinitialization [0.95, 1.0] (and after 1 step)
parameters after nsteps_neq steps [0.0, 0.6666666666666665]

I'm using a symmetric protocol for the alchemical_functions, so after nsteps_neq the electrostatics and sterics should both be 1.0 which they're not.

I found a strange fix to this behavior, which was to get an integrator variable (ex. integrator.getGlobalVariableByName('lambda') after every step taken. However strange thing was that if this was done immediately after reinitialization I got a segmentation fault.

#after adding the above fix
parameters ['lambda_electrostatics', 'lambda_sterics']
parameters before reinitialization [0.0, 0.6666666666666667]
parameters after reinitialization [0.0, 0.6333333333333334] (and after 1 step)
parameters after nsteps_neq steps [0.9999999999999998, 1.0]

It's not clear to me whether this is an integrator bug, an OpenMM bug, or intended behavior from OpenMM I just wasn't expecting. Since you folks are more familiar with the AlchemicalNonequilibriumLangevinIntegrator I thought I'd ask here first for thoughts on this.

I've incuded example scripts to replicate the above issue in param_issue.zip, which contains:

issue.py, reproduces the issue with the lambda_sterics and lambda_electrostatics being reset fix.py, is essentially the same code, with the integrator.getGlobalVariableByName('lambda') fix above to resolve the problem seg_fault.py, the same code again, but with a integrator.getGlobalVariableByName('lambda') line immediately after reinitializing.

[jchodera]

However when I coded this up I got behavior I didn't expect, namely that after reinitializing the lambda_sterics and lambda_electrostatics parameters were reset to 1.0.

This is expected OpenMM behavior. Whenever you reinitialize a Context, it should go back to the default global parameters set in the Forces, as if you had created a new Context.

I bet you could capture the global parameters with context.getParameters and then iterate over the map to context.setParameter after reinitializing.

[jchodera]

seg_fault.py, the same code again, but with a integrator.getGlobalVariableByName('lambda') line immediately after reinitializing.

This probably indicates a bug, and should probably be submitted as a bug report to OpenMM.

[jchodera]

I'm trying to reproduce this issue by running issue.py, since I'm concerned about whether the integrator global variables are reset when reinitializing the Context. I'm wondering if the output of my issue.py matches yours, however, since this doesn't match any of the output you pasted above:

arameters ['lambda_electrostatics', 'lambda_sterics']
parameters before reinitialization [0.0, 0.6666666666666667]
parameters after reinitialization [0.0, 0.6666666666666667] (and after 1 step)
parameters after nstepsNC steps [0.0, 0.6666666666666665]

[sgill2]

This is expected OpenMM behavior. Whenever you reinitialize a Context, it should go back to the default global parameters set in the Forces, as if you had created a new Context.

That makes sense. Thanks!

I'm trying to reproduce this issue by running issue.py, since I'm concerned about whether the integrator global variables are reset when reinitializing the Context. I'm wondering if the output of my issue.py matches yours, however, since this doesn't match any of the output you pasted above

@jchodera I must have copied some old output from a previous iteration of the script, I see that same output I get rerunning issue.py. In the way the code's written though, I still think, as you suggested, that there's a potential issue with the integrator variables being affected by reinitializing. The steric and electrostatics parameters should be updated each step based on the integrator's lambda variable, which isn't happening in the output above.

Does that referenced PR help address the behavior I'm seeing? Reading it, it wasn't clear to me what the surprising behavior mentioned was.