Error with restarting LAMMPS simulation with PLUMED and SHAKE

[ppegolo]

Dear developers,

there seems to be an issue with restarting a LAMMPS simulation that uses both PLUMED and the SHAKE algorithm.
When I try to restart a biased simulation from a LAMMPS restart file, I immediately receive the error Non-numeric box dimensions - simulation unstable. If the simulation had continued from the same point without restarting, I would not have encountered this error at this stage.

I have attached a ZIP archive with the LAMMPS (v. 23 Jun 2022) and PLUMED (v. 2.8.1) input files for reproducing the error.
The error appears when running step 4.
Does anybody have any idea if this is a bug of some sort?

Thanks,
Paolo

error_restart.zip

[invemichele]

There seems to be an issue with restart when a bias acts on the energy and the simulation is NPT.
Here test_restart.zip is a minimal setup to reproduce the error. The plumed input is:

ene: ENERGY
metad: METAD ARG=ene PACE=200 SIGMA=20 HEIGHT=1

Upon restart I get a huge spike in the pressure (from 380 to -10^227) and the simulation breaks.
There are no problems if a CV different from the energy is biased.

If I turn off the barostat and do NVT, I still get the pressure spike, but the simulation continues without issues.
In both cases, NPT and NVT, the kinetic energy after restart is identical, while the potential energy changes only slightly

[PabloPiaggi]

Hi @invemichele, we had never seen this error before. I believe it should be due to recent changes in either LAMMPS or Plumed. Did you test the same versions @ppegolo used? That is to say, LAMMPS (v. 23 Jun 2022) and PLUMED (v. 2.8.1). Thanks.

[invemichele]

yes, I used the latest LAMMPS stable (v. 23 Jun 2022) and let it download PLUMED (v. 2.8.1)

[ppegolo]

I also tried LAMMPS v. 29 Sep 2021 with PLUMED 2.7.4 with Michele's minimal setup and I get the same error.

[invemichele]

It might be connected to how the virial is handled in a restart...

Maybe related to shake? see here.
We should test if the jump in pressure happens also without the shake fix (probably smaller time step is needed with flexible water)

[ppegolo]

Yes, the error is related to shake. Removing shake (and reducing the time-step) makes the error disappear.

[invemichele]

I did the following test: add fix_modify shake1 virial no after fix shake1 and run with NVT (in NPT volume explodes)
With this setup, at restart there is no jump in the pressure.

[ppegolo]

Ah, I see. I'm not sure it's related, but in NVT simulations the Non-numeric box dimensions error does not show up at restart even without adding fix_modify shake1