Energy bug in PlumedForce

[stefdoerr]

Getting PLUMED energies a second time in the same step gives wrong results.

context = simulation.context
for i in range(10):
    simulation.step(100)
    ene1 = context.getState(getEnergy=True, groups={21}).getPotentialEnergy()
    ene2 = context.getState(getEnergy=True, groups={21}).getPotentialEnergy()
    assert ene1 == ene2, "Energy has changed!!!"

0.0 kJ/mol
10.0 kJ/mol
Traceback (most recent call last):
  File "/home/sdoerr/Work/pyacemd/debug_plumed/test.py", line 63, in <module>
    assert ene1 == ene2, "Energy has changed!!!"
AssertionError: Energy has changed!!!

Only the first energy matches the one reported by PLUMED itself if you dump out the bias of the force directly.

I attach a self-contained example which you can run: plumed_energy_bug.zip

This complicates things a lot if you add an energy reporter to the simulation, which is how it happened in my case, since I was requesting the energy once to get the total potential energy of the system and another time to get the PLUMED energy only and was thus getting wrong energies back.

[stefdoerr]

This seems to not happen with RESTRAINT PLUMED forces.

This here will wrongly change energies the second time you retrieve them:

phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
mtd: METAD ARG=phi,psi PACE=100 HEIGHT=10 SIGMA=0.1,0.1 BIASFACTOR=5 FILE=mtd.txt
PRINT ARG=mtd.bias STRIDE=100 FILE=bias.txt
FLUSH STRIDE=100

While this one here works fine:

pos: POSITION ATOM=7
rest: RESTRAINT ARG=pos.x,pos.y,pos.z KAPPA=100,100,100 AT=0,0,0
PRINT ARG=rest.bias FILE=bias.txt
FLUSH STRIDE=1

[stefdoerr]

I forgot to mention that to test this with my python script you need to compile #90 since it requires the temperature to be set.
But I believe it's not related to python and can probably be reproduced through the C++ layer as well.

[peastman]

What version are you using? This sounds like the bug that was fixed in #82. The fix is in the current release (2.0.1).

[tonigi]

Possibly related? #64

[stefdoerr]

I'm using the latest #90 build with the python fixes

…

On Wed, Nov 13, 2024, 19:50 Toni G ***@***.***> wrote: Possibly related? #64 <#64> — Reply to this email directly, view it on GitHub <#91 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB4RLAVQJUTVXG5LZBMJBKL2AOGPJAVCNFSM6AAAAABRWQCAP6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINZUGMZTGNZRGE> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[peastman]

I think this is inherent in your PLUMED script? You tell it to do metadynamics. When it evaluates the energy, it updates the bias so that future evaluations will produce different results.

[stefdoerr]

#64 is definitely the same issue, thanks Toni.
As the user mentions

an energy evaluation should not be counted as a step, right?

If what you are saying Peter is correct, if I never called getState it would not do metadynamics at all then?

[stefdoerr]

I found a (probably old) implementation of PLUMED in ASE https://wiki.fysik.dtu.dk/ase/_modules/ase/calculators/plumed.html
Generally it seems to do the same steps to get the energy and forces but there is a significant difference in that ASE caches the energy/force results inside of its Calculators if the conformation has not changed. So no matter how many times you call the energy/force evaluation, the energy/force calculation will only be done once.
https://github.com/hainm/ase/blob/master/ase/calculators/calculator.py#L458-L461

[tonigi]

Thank you for the work you are doing on this interface. <3

[peastman]

This seems to be an impedance mismatch. Here is how OpenMM expects integration to work.

1. Compute forces.
2. Update the positions and velocities.
3. Increment the step index.

PlumedForce is a Force. Everything happens inside PlumedForceImpl::computeForce(), which is expected to compute forces but not modify any state.

But PLUMED does update state. Here is the relevant code from the plugin.

openmm-plumed/openmmapi/src/PlumedForceImpl.cpp

Lines 143 to 148 in 95bfd46

	plumed_cmd(plumedmain, "prepareCalc", NULL);
	plumed_cmd(plumedmain, "performCalcNoUpdate", NULL);
	if (step != lastStepIndex) {
	plumed_cmd(plumedmain, "update", NULL);
	lastStepIndex = step;
	}

We call performCalcNoUpdate, which should not update the biases. Then we call update only if the step index has changed. (See https://www.plumed.org/doc-v2.9/developer-doc/html/_how_to_plumed_your_m_d.html.) But the update is happening at the wrong point in the step, as far as OpenMM is concerned. Incrementing the step index happens later in the step, after force computations are complete. So integrating a step doesn't immediately update the biases. Instead it causes them to be updated on the next force/energy computation, whether or not that is part of an integration step.

The correct solution, as far as OpenMM is concerned, is for PlumedForceImpl to implement updateContextState() and do the bias update there. That is typically at the start of the step, before computing forces. But that conflicts with PLUMED, which expects commands to be called in a particular order (first do the calculation, then perform the update).

[stefdoerr]

Interesting! Thanks so much for looking into it! I am curious though.

I imagine the computeForce is calculated once when I call the simulation.step(1) before the step is incremented as you said. Then I call computeForce manually a second time in my loop so now it gives the energy of the correct step. But when I call it a third time why does that modify the state again? It should not be able to enter that if clause and update PLUMED, correct?

for i in range(10):
    simulation.step(1) # Implicit computeForce call. Wrong energy
    ene1 = context.getState(getEnergy=True, groups={21}).getPotentialEnergy() # Correct energy
    ene2 = context.getState(getEnergy=True, groups={21}).getPotentialEnergy() # Wrong energy

[peastman]

None of the energies are wrong. They're all the correct energy at various points. The unexpected thing is when PLUMED updates the biases.

simulation.step(1)

This computes the forces, updates the positions based on them, and increments the step index.

ene1 = context.getState(getEnergy=True, groups={21}).getPotentialEnergy()

This computes the energy and returns it. In addition, it sees that the step index has changed since the last call, so it calls PLUMED with update, which causes it modify the biases. That happens after the energy has already been computed. PLUMED requires them to be called in that order, since the update needs the values of the CVs.

ene2 = context.getState(getEnergy=True, groups={21}).getPotentialEnergy()

This computes the energy again and returns it. The energy has changed, because PLUMED modified the biases after the last calculation. The step index has not changed, so it does not call update.